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FINAL  REPORT  ON  THE  SPECIAL  YEAR  IN  NUMERICAL  ANALYSIS,  1980-81 

1.  Introduction  //V /" (,'j/^  ~j^()  C / 

Each  year  the  Department  of  Mathematics  at  the  University  of  Maryland 
has  a  Special  Year  in  some  branch  of  mathematics.  This  past  year  the 
topic  was  Numerical  Analysis.  The  Special  Year  in  1980-81  was  jointly 
sponsored  by  the  Air  Force  Office  of  Scientific  Research  and  the  Mathematics 
Department.  The  main  goal  of  the  year  was  to  advance  the  state  of  the 
art  in  numerical  analysis  by  bringing  together  the  leading  experts  in  the 
field  for  formal  lectures  and  informal  discussions  of  recent  progress, 
current  problems,  and  future  trends.  We  placed  special  emphasis  on 
numerical  solution  of  partial  differential  equations,  global  continuation 
methods,  numerical  methods  in  statistics,  numerical  linear  algebra,  and 
numerical  problems  connected  with  special  functions.  The  activities  in 
statistics,  numerical  linear  algebra,  and  special  functions  took  place 
in  the  Fall  Semester  and  those  in  the  numerical  solution  of  partial  dif¬ 
ferential  equations  took  place  in  the  Spring.  In  the  area  of  global 
continuation  methods,  the  activities  were  spread  throughout  the  year.  A 
subsection  of  this  report  is  devoted  to  each  of  these  major  activities. 

Through  the  lectures  and  the  extensive  opportunities  for  informal 
discussions,  the  Special  Year  provided  an  excellent  opportunity  for  exchange 
of  information  and  ideas  between  the  members  of  the  large  and  active 
numerical  analysis  group  at  the  University  and  the  visiting  mathematicians. 

U'e  believe  the  benefits  will  be  substantial  both  to  the  University  and  vo 
the  international  numerical  analysis  community. 
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!.  Activities 

a.  Numerical  Solution  of  Partial  Differential  Equations 
A  wide  variety  of  important  problems  in  science  and  engineering  are 
formulated  as  initial-boundary  value  problems  for  partial  differential 
equations,  the  numerical  solution  of  which  is  one  of  the  most  important 
areas  of  numerical  analysis.  In  order  to  survey  the  central  problems  and 
trends  in  this  area,  we  invited  a  total  of  30  distinguished  visitors  to  the 
campus  during  the  Spring  Semester.  These  30  visitors  represented  nearly 
all  important  subfields  of  the  area.  The  lectures  by  the  visitors  were 
attended  by  Maryland  students  and  faculty  who  are  working  in  the  area 
or  have  an  interest  in  it,  members  of  the  numerical  analysis  community 
of  the  Washington  area  (e.g.,  from  the  Naval  Surface  Weapons  Center), 
and  of  course  the  other  visitors  in  residence  at  the  time.  In  addition 
to  the  lectures  there  was  ample  opportunity  for  informal  discussions. 

These  discussions  were  especially  fruitful  and  a  number  of  joint 
research  projects  have  grown  out  of  them.  Many  of  the  visitors  submitted 
written  versions  of  their  lectures.  These  range  from  extended  abstracts, 
to  systematic  survey  papers,  to  standard  research  papers.  This  collection 
of  papers  has  been  published  as  part  of  the  Lecture  Note  Series  of  the 
Mathematics  Department.  The  activities  in  this  area  were  loosely  divided 
between  finite  element  and  finite  difference  methods.  The  program  in 
Numerical  Solution  of  Partial  Differential  Equations  was  directed  by 
Professors  I.  Babuska,  T.-P.  Liu,  and  J.  Osborn. 

See  Attachment  A  for  a  list  of  participants  and  lectures. 
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b.  Numerical  Methods  in  Statistics,  Numerical  Linear  Algebra,  and 
Computation  of  Special  Functions 

The  fields  of  statistics,  linear  algebra,  and  special  functions  are 
usually  regarded  as  distinct  by  the  mathematical  community.  In  applications, 
however,  there  is  considerable  overlap:  each  field  leads  to  numerical 
and  other  problems  in  at  least  one  of  the  other  two.  The  increasing  speed 
and  capacity  of  modern  computers  is  bringing  more  and  more  of  these  problems 
into  the  realm  of  feasibility. 

In  order  to  generate  as  many  contacts  as  possible  between  research 
workers  in  the  three  fields,  it  was  decided  to  concentrate  activities 
into  a  single  conference.  This  was  held  from  October  2  to  October  8,  1980 
at  the  Adult  Education  Center  on  the  College  Park  campus.  The  conference 
was  advertised  in  the  Notices  of  the  American  Mathematical  Society,  SIAM 
News,  Bulletin  of  the  Institute  of  Mathematical  Statistics,  AMSTAT 
News,  and  the  Washington  Statistical  Society.  In  addition,  announcements 
of  the  conference  were  sent  directly  by  mail  to  a  list  of  approximately 
1650  individuals,  including  all  members  of  S.I.A.M.  resident  in  North  America. 
The  organizing  committee  consisted  of  F.  W.  J.  Olver,  G.  Stewart  and  G.  Yang. 

In  all,  132  individuals  registered  for  and  attended  the  conference. 

This  included  67  from  the  University  of  Maryland,  29  from  other  universities, 

18  from  government  agencies,  and  14  from  industry.  (Four  gave  no  affiliation.) 
The  attendees  from  the  University  of  Maryland  represented  11  different 
departments  or  programs. 


The  program  consisted  of  10  invited  1-hour  lectures  and  21  30-minute 
lectures.  Speakers  in  the  second  category  were  selected  by  the  organizing 
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committee  from  extended  abstracts  submitted  in  advance.  There  were  also 
sessions  for  15-minute  contributed  papers. 

See  Attachment  B  for  the  conference  program,  abstracts  of  the  1-hour 
invited  lectures,  abstracts  of  the  30-minute  selected  talks,  and  a  list  of 
registrants . 

c.  Homotopy  Continuation  Methods 

Homotopy  continuation  methods  are  directed  towards  solving  systems 
of  equations  in  situations  where  approximate  solutions  are  not  available 
and  quasi-newton  type  methods  fail.  Smooth  continuation  methods  were 
emphasized,  although  simplicial  methods  were  represented  by  Eaves  and 
Peitgen.  Keller's  and  Watson's  lectures  paid  special  attention  to 
applications  to  physical  problems.  Harrison  and  Smale  gave  lectures 
aimed  at  topics  relevant  to  theorical  understanding  of  families  of 
periodic  orbits  (Harrison)  and  the  number  of  steps  needed  to  implement 
methods  (Smale).  Interactions  with  visitors  have  led  to  two  papers  being 
written.  The  homotopy  program  was  organized  by  Professor  J.  A.  Yorke. 

See  Attachment  C  for  a  list  of  participants  and  lectures. 


LIST  OF  ATTACHMENTS 


ATTACHMENT  A:  List  of  Participants  and  Lectures  in  Numerical  PDE  Portion 
of  Special  Year.- 


ATTACHMENT  B:  Program,  Abstracts  of  the  1-Hour  Invited  Lectures,  Abstracts 
of  the  30-Minute  Selected  Talks,  and  List  of  Registrants 
for  the  Conference  on  Applications  of  Numerical  Analysis 
and  Special  Functions  in  Statistics,; 


ATTACHMENT  C:  List  of  Participants  and  Lectures  in  Global  Continuation 
Methods  Portion  of  Special  Year. 


ATTACHMENT  A:  List  of  Participants  and  Lectures  in  Numerical 


PDE  Portion  of  Special  Year 
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List  of  Participants 


PROFESSOR  GARTH  BAKER 
Harvard  University 

PROFESSOR  GARRETT  BIRKHOFF 
Harvard  University 

PROFESSOR  JAMES  BRAMBLE 
Cornell  University 

PROFESSOR  F.  BREZZI 
Universite  di  Pavia 
ITALY 

PROFESSOR  P.G .  CIARLET 
Universite  Pierre  et  Marie 
Curie,  FRANCE 

PROFESSOR  JIM  DOUGLIS,  JR. 
University  of  Chicago 

PROFESSOR  T.  DUPONT 
University  of  Chicago 

PROFESSOR  B.  ENQUIST 
University  of  California, 
Los  Angeles 

DR.  RICHARD  EWING 
Mobil  Field  Research 
Laboratory 

PROFESSOR  R.  FALK 
Rutgers  University 

PROFESSOR  P.  GARABEDIAN 
Courant  Institute 

PROFESSOR  J.  GLIMM 
Rockefeller  University 

PROFESSOR  AMIRAM  HARTEN 
Tel-Aviv  University 
ISRAEL 

PROFESSOR  LING  HSAIO 
Brown  University  and 
Academia  Sinica,  Peking 
PEOPLES  REPUBLIC  OF  CHINA 

PROFESSOR  P.  LAX 
Courant  Institute 


PROFESSOR  MITCHELL  LUSKIN 
Courant  Institute 

PROFESSOR  A.  MAJADA 
University  of  California, 
Berkeley 

PROFESSOR  J.  NITSCHE 
Inst.  Fur  Angewandte  Math. 
GERMANY 

PROFESSOR  J.T.  ODEN 
University  of  Texas, 

Austin 

PROFESSOR  J.  OLIGER 
Stanford  University 

PROFESSOR  A.  SCHATZ 
Cornell  University 

PROFESSOR  RIDGWAY  SCOTT 
University  of  Wisconsin 
Mathematics  Research  Center 

PROFESSOR  G.  STRANG 
Massachusetts  Institute 
of  Technology 

PROFESSOR  ROGER  TEMAM 
Universite  de  Paris 
FRANCE 

PROFESSOR  VIDAR  THOMEE 
Chalmers  University  of 
Technology 

PROFESSOR  LARS  WAHLBIN 
Cornell  University 

PROFESSOR  W.  WENDLAND 
Technische  Hochschule 
Darmstadt,  GERMANY 

PROFESSOR  B.  WENDROFF 
Los  Alamos  Scientific 
Laboratory 

PROFESSOR  MARY  WHEELER 
Rice  University 
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PROFESSOR  MILOS  ZLAMAL 
Technical  University 
CZECHLOSLOVAKIA 
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List  of  Lectures 


DATE 

SPEAKER 

Jan  16 

V. 

Thome e 

Jan  20 

V. 

Thome e 

Jan  20 

J. 

Bramble 

Jan 

22 

G. 

Birkhoff 

Jan 

27 

V. 

Thome e 

Jan 

29 

R. 

Scott 

Feb 

5 

M. 

Luskin 

Feb 

10 

A. 

Harten 

Feb 

10 

M. 

Luskin 

Feb 

10 

P. 

Lax 

Feb 

12 

R. 

Temam 

Feb 

12 

P. 

Ciarlet 

Feb 

17 

P. 

Ciarlet 

Feb 

17 

M. 

Wheeler 

Feb 

18 

R. 

Temam 

Feb 

19 

J. 

Nitsche 

Feb 

24 

J. 

Douglas 

Feb 

24 

J. 

Oden 

Feb 

26 

G. 

Baker 

TOPIC 

Single  step  methods  for  linear  differential  equations 
in  Banach  spaces,  PART  I 

Single  step  methods  for  linear  differential  equations 
in  Banach  spaces,  PART  II 

Remarks  on  Lagrange  multiplier  techniques  in  conjunc¬ 
tion  with  finite  element  approximations  in  various 
elliptic  problems 

Adapting  Courant-Friedrichs  Levy  to  the  1980' s 

Single  step  methods  for  linear  differential  equations 
in  Banach  spaces,  PART  III 

A  comparison  of  laboratory  experiments  with  a  model 
equation  for  water  waves 

Analysis  of  a  fractional  step  method  for  fluid  flow 
in  a  pipe 

On  random  choice  methods  for  hyperbolic  conservation 
laws 

On  a  finite  element  method  to  solve  the  critcality 
eigenvalue  problem  for  the  transport  equation 

Convergence  almost  everywhere  of  random  choice  schemes 

Variational  problems  in  mechanics  (plasticity)  PART  I 

Questions  of  existence  in  non  linear  elasticity 

Justification  of  the  von  Karman  equations 

Mixed  methods  for  miscible  displacement  problems 

Variational  problems  in  mechanics  (plasticity)  PART  II 

The  method  of  straightening  the  free  boundary  in 
moving  boundary  problems 

Numerical  simulation  of  flow  in  porus  media 

Analysis  of  some  contact  problems  in  nonlinear 
elasticity 

Spectral  approximation  in  Riemannian  geometry 
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List  of  Lectures  (-continued) 


DATE 

SPEAKER 

TOPIC 

Feb 

26 

R. 

Ewing 

Alternating  directional  multistep  procedure  for 
nonlinear  parabolic  P.D.E.'s 

Feb 

27 

J. 

Nitsche 

A  remarkable  approximation  property  of  finite  ele¬ 
ments  and  its  consequences 

Mar 

3 

B. 

Enquist 

Radiation  boundary  conditions  at  computational 
boundaries 

Mar 

3 

B. 

Wendrof f 

Defect  corrections,  multigrids  and  selected  appli¬ 
cations,  PART  I 

Mar 

5 

L. 

Wahlbin 

On  maximum  norm  estimates  in  finite  element  methods 
PART  I 

Mar 

6 

W. 

Wendland 

Asymptotic  convergence  of  boundary  element  methods 

Mar 

10 

B. 

Enquist 

Flux  splittings  in  compressible  flow  computations 

Mar 

10 

L. 

Wahlbin 

On  maximum  norm  estimates  in  finite  element  methods 
PART  II 

Mar 

12 

B. 

Wendroff 

Defect  corrections,  multigrids  and  selected  appli¬ 
cations,  PART  II 

Mar 

12 

W. 

Wendland 

Integral  equation  methods  for  mixed  boundary  value 
problems 

Mar 

12 

L. 

Wahlbin 

On  maximum  norm  estimates  in  finite  element  methods 
PART  III 

Mar 

24 

A. 

Majda 

Vortex  methods  in  fluid  flow 

Mar 

24 

T. 

Dupont 

Mesh  modification  in  finite  element  methods 

Mar 

26 

J. 

Glimm 

Hydrodynamics  without  diffusion:  Theory,  computation 
and  application,  PART  I 

Mar 

31 

A. 

Schatz 

Singular  functions  in  the  finite  element  method 

Apr 

1 

A. 

Majda 

A  theory  for  Mach  Stern  formation  in  reacting  shock 
fronts 

Apr 

2 

F. 

Brezzi 

Finite  dimensional  approximation  of  nonlinear  prob¬ 
lems,  PART  I 

Apr 

7 

A. 

Schatz 

Boundedness  in  L  of  the  Rite  projection 

<X> 

Apr  7 


M.  Zlamal 


Boundedness  in  L  of  the  Rite  prelection 

CX> 

Galerkin-f inite  element  met’ncas  for  the  solution  of 
nonlinear  evolution  equations,  PART  I 
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List  of  Lectures  (continued) 


DATE 

SPEAKER 

TOPIC 

Apr 

8 

J. 

Glimm 

Hydrodynamics  without  diffusion:  Theory,  computation 
and  application,  PART  II 

Apr 

9 

L. 

Hsiao 

Overtaking  of  shock  waves  in  steady  two  dimensional 
supersonic  flows 

Apr 

9 

J. 

Oliger 

Adaptive  difference  methods  for  time  dependent  prob¬ 
lems 

Apr 

10 

J. 

Glimm 

Mathematical  aspects  of  quantum  field  theory 

Apr 

14 

J. 

G 1  imm 

Hydrodynamics  without  diffusion:  Theory,  computet icr. 
and  application,  PART  III 

Apr 

14 

M. 

Zlamal 

Galerkin-f inite  element  metnods  for  the  solution  of 
nonlinear  evolution  equations,  PART  II 

Apr 

15 

R. 

Falk 

A  mixed  finite  element  method  for  the  simply  sup¬ 
ported  plate  problem 

Apr 

16 

G. 

Strang 

Optimal  design 

Apr 

16 

F. 

Brezzi 

Finite  dimensional  approximation  of  nonlinear  prob¬ 
lems,  PART  II 

Apr 

21 

F. 

Brezzi 

Finite  dimensional  approximation  of  nonlinear  prob¬ 
lems,  PART  III 

Apr 

21 

P. 

Garabedian 

Numerical  analysis  of  equilibria  with  islands  in 
magnet  ohyd  rodynamics 

Apr 

23 

M. 

Zlamal 

Galerkin-f inite  element  methods  for  the  solution  of 
nonlinear  evolution  equations,  PART  III 

Apr 

24 

F. 

Brezzi 

Finite  dimensional  approximation  of  nonlinear  prob¬ 
lems,  PART  IV 

Apr 

28 

M. 

Zlamal 

Galerkin-f inite  element  methods  for  the  solution  of 
nonlinear  evolution  equations,  PART  IV 

ATTACHMENT  B:  Program,  Abstracts  of  the  1-Hour  Invited  Lectures, 
Abstracts  of  the  30-Minute  Selected  Talks,  and  List 
of  Registrants  for  the  Conference  on  Applications  of 
Numerical  Analysis  and  Special  Functions  in  Statisti 


program  for 

THE  CONFERENCE  ON  APPLICATIONS  OF  NUMERICAL 
ANALYSIS  AND  SPECIAL  FUNCTIONS  IN  STATISTICS 


held  at 


Adult  Education  Center 
UNIVERSITY  OF  MARYLAND 
College  Park,  Maryland 

October  2  -  8 ,  1980 


sponsored  by 


U.S.  Air  Force  Office  of  Sc i ent i f i c  Resea rc h 
Department -of  Mathematics,  University  of  Maryland 


organ  t 


i  n  g 


committee 


Grace  Yang  The  Mathematical  Statistics  Program 

Frank  Olver  Institute  for  Physical  Science  and  Technology 

G  .  VI .  Stewart  Department  of  Computer  Science 


l-i 


special 


year 


numerical 


analysis: 


i 


n 


This  conference  is  part  of  a  Special  Year  in  Numerical  Analysis  sponsored 
by  the  U.S.  Air  Force  Office  of  Sc  i  ent  i  f  ic  Research  and  the  Department  of 
Mathematics,  University  of  Maryland.  Special  Year  visitors  are  listed  be¬ 
low. 

For  further  information  contact  Professors:  P.  Wolfe 

Y.  Yorke 
I .  Babuska 


SPECIAL  YEAR  VISITORS: 


R.  Askey* 

Spring  G.  Birkhoff 

A.  Harten 

G.  Strang 

W.  Coo/ 

J .  Bramb I e 

H.  Keller 

R.  Temam 

J.  Dennis* 

F.  Brezzi 

H.-O.  Kreiss 

V.  Thomee 

W.  Gautschi* 

A.  Chorin 

P.  Lax 

W.  Wendlund 

J.  Harrison 

P.  Ciarlet 

M.  Luskin 

B.  Wendroff 

S.  Karl i n* 

J.  Douglas,  Jr. 

A.  Majda 

M.  Wheeler 

P.  Krishnaiah* 

T.  Dupont 

J.  Nitsche 

M .  Z 1 ama 1 

J .  Lyness* 

B.  Enquist 

J .  Oden 

C.  Paige* 

R.  Ewing 

J .  01 iger 

H.O.  Peitgen 

R.  Falk 

A.  Schatz 

J .  Snel 1* 

P.  Garabedian 

R.  Scott 

L.  Watson 

J .  G 1 i mm 

S .  Sma 1 e 

^Conference  speakers. 


The  registration  booth  for  the  Conference  is  located  on  the  Main  Concourse 
in  the  Adult  Education  Center  and  will  be  manned  as  follows: 

Wednesday  October  I  7:30  -  9:30  P.M. 

Thursday  October  2  8:30  -  12:00  A.M. 

Friday  October  3  8:30  -  12:00  A.M. 

Saturday  October  4  8:30  -  12:00  A.M. 

Sunday  October  5  7:30  -  9:30  P.M. 

Monday  October  6  8:30  -  12:00  A.M. 

Tuesday  October  7  8:30  -  12:00  A.M. 

Message  Board 

There  will  be  a  message  board  located  at  the  Conference  Room,  #1105,  for 
your  convenience. 


Additional  Information 


The  names  of  the  morning  and  afternoon  chairpersons  and  the  names  and 
titles  of  the  speakers  giving  contributed  papers  will  be  published  sep¬ 
arately,  since  they  are  not  available  at  the  time  of  this  printing. 

There  will  also  be  available  a  phamplet  containing  information  on  local 
transportation  and  sightseeing,  and  restaurants. 


Thursday  10-2-80  Rm  1105 


9:00  -  9:30  OPENING  OF  THE  CONFERENCE:  W.E.  KIRWAN 

(Chairman,  Dept,  of  Mathematics) 

WELCOMING  ADDRESS:  FRANK  J .  KERR 

(Provost,  Division  of  Mathematical 
&  Physical  Sciences  &  Engineering ) 

MORNING  SESSION: 

9:30  -  10:30  C.C.  PAIGE,  McGill  University,  CANADA 

"The  General  Gauss-Markov  Model  and  the  Singular  Value 
Decomposition" 

10:30  -  11:00  Coffee 

11:00-  11:30  GEORGE  CYBENKO,  Tufts  University 

"The  Efficient  Solution  by  Orthogonal izat ion  of  Linear 
Prediction  Problems  for  Stationary  Time  Series" 

11:30  -  12:00  JAMES  A.  CADZOW,  Virginia  Polytechnic  Inst. 

"Autoregressive-Moving  Average  Spectral  Estimation:  A 
New  Effective  Model ing  Procedure" 

AFTERNOON  SESSION: 

1:30-2:30  S.  KARLIN,  Stanford  University 

"A  Diffusion  Stochastic  Model  of  Mathematical  Genetics 
Involving  Airy  Functions" 

2:30  -  3:00  Coffee 

3:00-3:30  S.K.  KATTI,  University  of  Missouri 

Topic:  Infinite  Divisibility 

3:30  -  4:00  K.O.  BOWMAN,  Union  Carbide  Corporation 

"Models  for  Approximating  the  Percentage  Points  of 
Di str i but  ions" 


4:00  -  4:30 


Contributed  Papers 


Friday  10-3-80  RM  1105 


MORNING  SESSION: 

9:00-10:00  J.E.  DENNIS,  Rice  University 

"Inside  Optimization  Routines" 

10:00  -  10:30  Coffee 

10:30  -  11:00  FRANKLIN  T.  LUK,  Cornell  University 

"The  Communal ity  Problem  for  Stieltjes  Matrices" 

11:00  -  11:30  ROBERT  B.  DAVIES,  University  of  California ,  Berkeley 
"Maximum  Likelihood  Estimation" 

11:30  -  12:00  Contributed  Papers 

AFTERNOON  SESSION: 

1:30  -  2:30  P.R.  KRISHNAIAH,  University  of  Pittsburgh 

"Computations  of  Multivariate  Distributions" 

2:30  -  3:00  Coffee 

3:00-3:30  LOUIS  KATES 

"The  Zonal  Polynomials  of  Multivariate  Analysis  as 
Special  Functions" 

3:30  -  4:00  KEVIN  W.J.  KADELL,  University  of  Wisconsin 

"The  Selberg  Distribution" 

4:00  -  4:30  Contributed  Papers 

Saturday  10-4-80  R’^  1123 

MORN  I NG  SESSION: 

9:00  -  10:00  J.N.  LYNESS,  Argonne  National  Laboratory 

"The  Calculations  of  Trigonometric  Fourier  Coefficients" 

10:00  -  10:30  Coffee 


(SATURDAY  SESSION  CONTINUED  ON  NEXT  PAGE) 


10:30  -  11:00  PAUL  SPECKMAN,  University  of  Oregon 

"Spline  Smoothing  and  Optimal  Rates  of  Convergence  in 
Nonparametr ic  Regression  Models" 

11:00  -  11:30  MICHAEL  GHIL,  Courant  Institute 

"A  Stochastic-Dynamic  Model  for  Global  Atmospheric  Mass 
Field  Statistics" 

Monday  10-6-80  RM  1105 

MORNING  SESSION: 

9:00  -  10:00  J.L.  SNELL,  Dartmouth  University 

"Random  Walks  and  Electric  Networks" 

10:00  -  10:30  Coffee 

10:30  -  11:00  ROY  S.  FREEDMAN,  Hazeltine  Corporation 

"Random  Walks  and  Statistical  Commun ication" 

11:00  -  11:30  DAVID  THOMSON,  Bell  Laboratories 

"Applications  of  Spheriodal  Wave  Functions  to  Time  Series 
Analysi s" 

11:30  -  12:00  Contributed  Papers 

AFTERNOON  SESSION: 

1:30  -  2:30  G.W .  STEWART,  University  of  Maryland 

"Matrix  Perturbation  Theory  and  Linear  Regression" 

2:30  -  3:00  Coffee 

3:00-3:30  E.J,  WEGMAN,  Office  of  Naval  Research 

"On  Computer  Architectures  for  Statistical  Algorithms" 

3:30  -  4:00  ASH  IS  SEN  GUPTA,  Stanford  University 

"On  the  Applications  of  Special  Function-.  :n  To.'  ■  i 
Standardized  Generalized  Variances  of  Mui i ivariai**  Nor¬ 
mal  Populations  of  Possibly  Diffcrcr. i  Dimensions” 


(MONDAY  SESSION  CONTINUED  ON  NEXT  PAGE) 


4:00  -  4:30  Contributed  Papers 


EVENING: 

7:30  -  9:30  Wine  and  cheese  garden  party  at  the  Rossborough  Inn 
Tuesday  10-7-80  RM  1105 
MORNING  SESSION: 

9:00  -  10:00  W.J,  CODY,  Argonne  National  Laboratory 

"Preliminary  Report  on  Software  for  the  Modified  Bessel 
Functions  of  the  First  Kind" 

10:00  -  10:30  Coffee 

10:30  -  11:00  N.M.  TEMME,  Mathematisch  Centrum ,  THE  NETHERLANDS 

"Incomplete  Gamma  Functions,  Numerical  and  Asymptotical 
Aspects  for  Evaluation  and  Inversion" 

11:00  -  11:30  RODERICK  WONG,  University  of  Manitoba ,  CANADA 

"Some  Applications  of  Asymptotics  in  Statistics" 

11:30  -  12:00  Contributed  Papers 

AFTERNOON  SESSION: 

1:30  -  2:30  R.A.  ASKEY,  University  of  Wisconsin 

"Gamma  and  Beta  Functions  From  Euler  to  Selberg  and  Their 
Orthogonal  Complements" 

2:30  -  3:00  Coffee 

3:00  -  3:30  CHARLES  F.  DUNKL,  University  of  Virginia 
"Discrete  Orthogonal  Polynomials" 

3:30  -  4:00  MARIETTA  J.  TRETTER,  The  Pennsylvania  State  University 

"Absolute  Error  Bounds  for  Edgeworth  Asymptotic  Expansion 

4:00  -  4:30  Contributed  Papers 


Wednesday  10-8-80  RM  1105 


MORNING  SESSION: 

9:00  -  10:00  W.  GAUTSCHI ,  Purdue  University 

"Special  Functions:  Computational  Considerations" 

10:00  -  10:30  Coffee 

10:30  -  11:00  DONALD  E.  AMOS,  Sandia  National  Laboratories 

"Computations  of  the  Central  and  Noncentral  F 
Di str i but i ons" 

11:00  -  11:30  WALTER  R.  NUNN,  Center  for  Naval  Analysis 
"The  Laguerre  Transform" 


PROGRAM  SUPPLEMENT: 
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THE  CONFERENCE  ON  APPLICATIONS  OF  NUMERICAL 
ANALYSIS  AND  SPECIAL  FUNCTIONS  IN  STATISTICS 

\ 

CHAIRPERSONS  FOR  THE  CONFERENCE: 


Thursday,  Oct.  2, 

1980 

Morn i ng : 
Afternoon : 

Professor 

Professor 

Diane  O'Leary 
Peter  Wolfe 

Friday, 

Oct.  3, 

1980 

Morn i ng : 
Afternoon : 

Professor 
Prof essor 

G.W.  Stewart 
Paul  Smith 

Saturday,  Oct.  4, 

1980 

Morning: 

Professor 

B.  Ke 1 logg 

Monday, 

Oct.  6, 

1980 

Morn i ng : 
Afternoon : 

Professor 

Professor 

G.  Yang 

S.  Kotz 

Tuesday, 

Oct.  7, 

1980 

Morn i ng : 
Afternoon : 

Professor 

Professor 

F.  Olver 

B.  Carlson 

Wednesday  Oct.  8, 

1980 

Morn i ng : 

Professor 

J .  Kei 1  son 

PARTIAL  LIST  OF  CONTRIBUTED  PAPERS: 
Dr.  E.  Cu+bi I  I 

David  Taylor  Naval  Ship  R  &  D  Center 


Mon.  Dr.  Alexander  S.  Elder 

IC-6-80  Aberdeen  Proving  Ground 

4:00 


"Ascending  and  Asymptotic  Series  for  Squares, 
Products  and  Cross  Products  of  the  Modified 
Bessel  Functions” 


Dr.  Jerry  Leon  Fields 
University  of  Alberta 


Mon.  Dr.  James  W.  Long  ley 

1 0-6-80 
11:45 

Fri.  Dr.  Clifford  Spiegelman 

10-5-80  National  Bureau  of  Standards 
!  1  : 30 

Mon.  Dr.  James  N.  Walbert 

10-6-30 

11:30 


Topic:  Convergence  of  an  explicit  sequence  of 
rational  approximations  to  the  hypergeometr i c 
functions 

.a 


F(b)  = 


q+ 


in  the  region 
D  =  {v  :  |arg  vj  <  tt. 


p  fal’*'-'aq+l  -T) 

> . % 


|arg(l+v)|  <  tt,  v  t  C-l,03}. 


"Modified  Gram-Schmidt  Process  Versus  Classical 
Gram  Schmidt" 


"An  Algorithm  for  Minimizing  an  Implicitly 
Restricted  Objective  Function" 

(with  Dr.  William  J.  Studden,  Purdue  Uniy.) 

"Use  of  a  Continued  Fraction  to  Evaluate  the 
Exponential  Integral  in  the  Complex  Plane" 


Thiurs. 

10-2-80 

4-4:30 


SPECIAL  HALF-HOUR  TALK: 
Dr.  Richard  Heilberger 


"The  Design  and  Construct  of  Test  Data  Sets  for 
Regression  Procedures"  (with  Dr.'s  Paul  F. 

Ve Neman  and  Agelia  Ypellar) 


PROGRAM 


SYNOPSIS 


PRINCIPAL  SPEAKERS 

(Fifty-Minute  Lectures) 


R. A.  ASKEY 

University  of  Wisconsin 

"Gamma  and  Beta  Functions  From  Euler  to  Selberg 
and  Their  Orthogonal  Polynomials" 

W.J.  CODY 

Argonne  National  Laboratory 

"Preliminary  Report  on  Software  for  the  Modified 
Bessel  Functions  of  the  First  Kind" 

J.E.  DENNIS 

Rice  University 

"Inside  Optimization  Routines" 

W.  GAUTSCH 1 

Purdue  University 

"Special  Functions:  Computational  Considerations 

S.  KARLIN 

Stanford  University 

"Various  Methods  for  Calculating  Family  Correla¬ 
tions  With  Variable  Family  Size" 

P.R.  KRISHNAIAH 

University  of  Pittsburgh 

"Computations  of  Multivariate  Distributions" 

J.N.  LYNESS 

Argonne  National  Laboratory 

"The  Calculation  of  Tr igunometr ic  Fourier  Coef¬ 
ficients" 

C.C.  PAIGE 

McGill  University 

"The  Genera!  Gauss-Markov  Model  and  the  Singular 
Value  Decomposition" 

J.L.  SNELL 

Dartmouth  University 

"Random  Walks  and  Electric  Networks" 

G.W.  STEWART 

University  of  Maryland 

"Matrix  Perturbation  Theory  and  Linear  Regression' 

DONALD  E.  AMOS 

Sandia  National  Laboratories 
K.O.  BOWMAN 

Union  Carbide  Corporation 

JAMES  A.  CADZOW 
Virginia  Polytechnic  Inst. 


INVITED  SPEAKERS 

(Ha I f-Hour  Talks) 

"Computation  of  the  Central  and  Noncentral  F 
Di str i but  ions" 

"Models  for  Approximating  the  Percentage  Points  of 
D i str i but  ions" 

"Autoregress i ve-Mov i ng  Average  Spectral  Estimation 
A  New  Effective  Modeling  Procedure" 
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\ 


GEORGE  CYBENKO 
Tufts  University 


ROBERT  B.  DAVIES 
Univ.  of  California,  Berkeley 

CHARLES  F.  DUNKL  "Discrete  Orthogonal  Polynomials" 


University  of  Virginia 

ROY  S.  FREEDMAN 

Hazeltine  Corporation 

"Random  Walks  and  Statistical  Communication" 

MICHAEL  GHIL 

Courant  Institute 

"A  Stochastic-Dynamic  Model  for  Global  Atmospheric 
Mass  Field  Statistics" 

ASH  IS  SEN  GUPTA 

Stanford  University 

"On  the  Applications  of  Special  Functions  in  Tests 
for  Standardized  Generalized  Variances  of  Multi¬ 
variate  Normal  Populations  of  Possibly  Different 
Dimensions" 

KEVIN  W.J.  KADELL 

University  of  Wisconsin 

"The  Selberg  Distribution" 

LOUIS  KATES 

"The  Zonal  Polynomials  of  Multivariate  Analysis  as 
Special  Functions" 

S.K.  KATTI 

University  of  Missouri 

Topic:  Infinite  Divisibility 

FRANKLIN  T.  LUK 

Cornell  University 

"The  Communal ity  Problem  for  Stieltjes  Matrices" 

WALTER  R.  NUNN 

Center  for  Naval  Analysis 

"The  Laguerre  Transform" 

PAUL  SPECKMAN 

University  of  Oregon 

"Spline  Smoothing  and  Optimal  Rates  of  Convergence 
in  Nonparametr ic  Regression  Models" 

N.M.  TEMME 

Mathematisoh  Centrum , 

THE  NETHERLANDS 

"Incomplete  Gamma  Functions,  Numerical  and  Asy-ptot 
cal  Aspects  for  Evaluation  and  Inversion" 

DAVID  THOMSON 

Bell  Laboratories 

"Applications  of  Spheroidal  Wave  Functions  tc  Time 
Series  Analysis" 

''ARIETTA  J.  TRETTER 

Tnc  Pennsylvania  State  Univ. 

"Absolute  Error  Bounds  for  Edgewcrtn  Asymp.ofi. 
Expansions" 

"The  Efficient  Solution  by  Orthogonal izat ion  of 
Linear  Prediction  Problems  for  Stationary  Time 
Series" 

"Maximum  Likelihood  Estimation" 
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E.J.  WEGMAN  "On  Computer  Architectures  for  Statistical 

Office  of  Naval  Research  Algorithms" 

RODERICK  WONG  "Some  Applications  of  Asymptotics  in  Statis- 

University  of  Manitoba  tics" 


TENTATIVE  SCHEDULE 


This  is  a  tentative  schedule  and  is  subject 

to  change. 

Thursday 

10-2-80 

9:00-9:30 

A.M. 

1 NTR0DUCT 1  ON 

9:30-10:30 

A.M. 

C.C.  Paige 

10:30-1 1 :00 

A.M. 

Coffee 

1 1 : 00- 12:00 

A.M. 

George  Cybenko 
James  A.  Cadzow 

1 : 30-2: 30 

P.M. 

S .  Karlin 

2:30-3:00 

P.M. 

Coffee 

3:00-4:00 

P.M. 

S.K.  Katti 

K.0.  Bowman 

4:00-4:30 

P.M. 

Contributed  Papers 

Fr i day 

10-3-80 

9:00-10:00 

A.M. 

J.E.  Dennis 

10:00-10:30 

A.M. 

Coffee 

10:30-1 1 : 30 

A.M. 

Frankl in  T.  Luk 
Robert  B.  Davies 

o 

1 

rsi 

o 

o 

A.M. 

Contributed  Papers 

1 : 30-2 : 30 

P.M. 

P.R.  Krishnaiah 

2:30-3:00 

P.M. 

Coffee 

3:00-4:00 

P.M. 

Louis  Kates 

Kevin  W. J .  Kadel 1 

4:00-4:30 

P.M. 

Contributed  Papers 

Saturday 

10-4-80 

9:00-10:00 

A.M. 

J.N.  Lyness 

10:00-10:30 

A.M. 

Coffee 

10:30-11:30 

A.M. 

Paul  Speckman 

Michael  Ghi I 


Monday 


J .L.  Snel I 


10-6-80  9:00-10:00  A.M. 


Tuesday 


Wednesday 


10:00-10:30  A.M. 

Coffee 

10:30-11 : 30  A.M. 

Roy  S.  Freedman 

David  Thomson 

1  1  -.30-12:00  A.M. 

Contributed  Papers 

1:30-2:30  P.M. 

G.W.  Stewart 

2:30-3:0 0  P.M. 

Coffee 

3:00-4:00  P.M. 

E.J.  Wegman 

Ash  is  Sen  Gupta 

4:00-4:30  P.M. 

Contributed  Papers 

t 0-7-80  9:00-10:00  A.M. 

W.J.  Cody 

10:00-10:30  A.M. 

Coffee 

10:30-1 1 :30  A.M. 

N.M.  Temme 

Roderick  Wong 

1  1 : 30- (2:00  A.M. 

Contributed  Papers 

1:30-2:30  P.M. 

R. A.  Askey 

2:30-3:00  P.M. 

Coffee 

3:00-4:00  P.M. 

Charles  F.  Dunkl 
Marietta  J.  Tretter 

4:00-4:30  P.M. 

Contributed  Papers 

10-8-80  9:00-10:00  A.M. 

W.  Gautschi 

10:00-10:30  A.M. 

Cof  fee 

10:30-1 1:30  A.M. 

Dona  Id  E .  Amos 
Walter  R.  Nunn 
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R.A.  Askey 

University  of  Wisconsin 

"Gamma  and  Beta  Functions  From  Euler  to  Selberg 
and  Their  Orthogonal  Polynomials" 

Many  of  the  classical  orthogonal  polynomials  first  arose  in  a  probabilistic 
or  statistical  setting.  Lagrange  discovered  Legendre  polynomials  and  their 
recurrence  relation  while  solving  a  discrete  probability  problem.  Laplace 
used  Hermite  polynomials  extensively  in  his  book  on  probability  theory. 
Fisher  rediscovered  the  discrete  Chebychev  polynomials  while  fitting  rain¬ 
fall  data.  Fisher's  representation  for  these  polynomials  was  different 
than  Chebychev 's  and  could  have  led  to  the  discovery  of  an  important  set  of 
orthogonal  polynomials  related  to  the  6-j  symbols  of  angular  momentum  theory 
if  anyone  had  seriously  looked  at  his  work.  From  these  orthogonal  polynom¬ 
ials  it  is  easy  to  find  the  three  term  recurrence  relation  for 


which  was  used  in  the  first  proof  of  the  irrationality  of  5(3).  In  a  com¬ 
pletely  different  field,  statistical  mechanics,  R.J.  Baxter  has  recently 
solved  another  two  dimensional  model  (the  hard  hexagon)  and  he  needed  the 
Rogers-Ramanu jan  identities  to  compute  a  limit  associated  with  phase  transi¬ 
tions.  These  identities  were  discovered  by  Rogers  while  studying  some 
polynomials  orthogonal  with  respect  to  measures  that  generalize  the  symmetric 
beta  and  normal  distributions.  A  brief  outline  of  these  beta  functions  will 
be  given  and  then  similar  integrals  in  several  variables  will  be  considered. 
After  work  of  Wishart,  Fisher,  Hsu,  Wilks  and  Ingham  in  statistics  and  Siegel 
in  number  theory,  the  first  real  break-through  was  made  by  A.  Selberg  in  1944, 
but  his  work  was  lost  for  almost  thirty-five  years.  Mehta  and  Dyson  extended 
Wishart’s  work  to  other  classes  of  matrices  and  came  up  with  a  beautiful  con¬ 
jectured  extension  of  the  normal  integral.  This  conjecture  is  easy  to  prove 
from  Selberg' s  integral.  Many  new  conjectures  have  been  formulated  in  the  las 
year  or  so.  A  few  of  these  will  be  mentioned. 
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W.J.  Cody* 

Argonne  National  Laboratory 

"Preliminary  Report  on  Software  for  the  Modified 
Bessel  Functions  of  the  First  Kind" 

In  our  experience  programs  which  evaluate  Bessel  functions  of  various  kinds 
are  requested  more  frequently  than  programs  for  any  other  special  functions. 
This  is  partially  because  of  the  importance  of  these  functions,  and  partially 
because  of  the  lack  of  high-quality  transportable  software  for  their  evalua¬ 
tion.  This  report  on  the  modified  Bessel  functions  of  the  first  kind  is  the 
first  of  a  series  of  projected  reports  surveying  available  Bessel  function 
software  and  laying  the  foundations  for  the  development  of  a  collection  of 
transportable  Bessel  function  programs. 

After  brief  discussions  of  relevant  analytic  properties  of  the  Bessel  func¬ 
tions,  important  computational  algorithms  derived  from  them,  and  desirable 
properties  of  good  numerical  software,  we  give  capsule  appraisals  of  eleven 
diverse  contemporary  programs  or  program  packages  for  the  I  Bessel  functions. 
We  then  describe  a  modification  of  one  of  the  more  promising  programs  to 
improve  its  performance  and  extend  its  capabilities.  Finally,  this  extended 
program  and  one  other  with  similar  capabilities  are  examined  in  greater 
detail  to  determine  whether  they  are  candidates  for  inclusion  in  the  proposed 
package. 


John  E.  Dennis 
Rice  University 

"Inside  Optimization  Routines" 


Applied  statisticians  often  find  library  subroutines  for  unconstrained  minimiza¬ 
tion  useful.  This  talk  will  attempt  to  explain  the  ideas  implemented  in  the 
best  routines.  We  will  also  mention  some  current  optimization  software  research 
directions . 


'"Work  performed  under  the  auspices  of  the  U.S.  Department  of  Energy. 
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Walter  Gautschl 
Purdue  University 

"Special  Functions:  Computational  Considerations" 

We  discuss  computational  aspects  of  infinite  series  and  continued  fractions 
in  the  context  of  evaluating  special  functions,  in  particular,  probability 
distribution  functions.  Questions  of  major  concern,  with  regard  to  infinite 
series,  are  internal  cancellation  of  terms  and  stopping  rules  for  truncat¬ 
ing  the  infinite  series.  We  introduce  appropriate  terminology  and  theory, 
and  give  applications  to  certain  power  series  related  to  the  incomplete 
gamma  function.  We  then  recall  briefly  how  various  types  of  continued  frac¬ 
tions  arise  through  correspondence  (or  association)  with  formal  power  series, 
or  via  second-order  linear  difference  or  differential  equations.  We  advocate 
Euler's  method  of  computing  a  continued  fraction  as  infinite  series.  A  new 
theorem  is  presented  concerning  the  convergence  behavior  of  continued  frac¬ 
tions  with  real  elements,  and  we  show  how  it  can  be  used  to  explain  the 
empirically  known  effectiveness  of  Legendre's  continued  fraction  for  the 
complementary  incomplete  gamma  function  T(a,x),  considering  that  conver¬ 
gence,  in  theory,  is  only  sublinear.  We  also  draw  attention  to  the  compu¬ 
tational  advantages  of  a  continued  fraction  of  Schlomilch  for  the  incomplete 
gamma  function  y(a,x). 


J.N.  Lyness 

Argonne  National  Laboratory 

"The  Calculation  of  Trigonometric  Fourier 
Coefficients" 

A  technique  for  the  numerical  approximation  of  sets  of  Trigonometric  Fourier 

•  f  1  2Trirx  « 

coefficients  I  f(x)e  dx;  r  =  0,1,2,...  based  on  a  common  set  function 

J  0 

values  f(x),  i  =  l,2,...,m  was  described.  The  underlying  theory  is  based 
on  subtracting  out  an  approximation  to  the  truncated  Euler  expansion  which  can 
be  integrated  analytically.  The  method  is  restricted  to  functions  having  a  high 
degree  of  continuity,  but  can  be  used  when  only  irregularly  spaced  function 
values  are  available. 

The  calculation  of  individual  Fourier  Coefficients  of  an  analytic  function  by 
using  contour  integration  in  the  complex  plane  was  also  discussed  briefly. 


C.C.  Paige 
McGill  University 
CANADA 


"The  General  Gauss-Markov  Model  and  the  Singular 
Value  Decomposition" 

The  problem  of  finding  the  best  linear  unbiased  estimate  §  (BLUE)  for  the 

2 

general  Gauss-Markov  linear  model  (y,XB,a  ft)  may  be  formulated  as  the  con¬ 
strained  linear  least  squares  problem 

T 

(1)  minimize  u  u  subject  to  y  =  Xg  +  Lu, 

T 

where  LL  =  ft  is  the  Cholesky  decomposition  of  the  given  nonnegative  definite 
symmetric  matrix  ft.  When  ft  is  positive  definite  L  is  nonsingular,  and  the 
singular  value  decomposition  of  L  ^X  could  in  theory  be  used  to  solve  this 
problem.  However,  such  an  approach  would  not  in  general  be  numerically  reliable, 
and  is  not  clearly  defined  when  L  is  singular  as  can  happen  in  practice. 

A  simultaneous  decomposition  of  L  and  X  is  suggested  which  is  based  on 
numerically  reliable  orthogonal  transformations  and  leads  immediately  to  the 
solution  of  (1).  This  decomposition  is  valid  for  all  L  and  X  with  the 
same  number  of  rows,  and  when  L  is  nonsingular  it  immediately  gives  the 
singular  value  decomposition  of  L  1X,  but  without  using  the  inverse  of  L. 

Thus  the  decomposition  is  an  appealing  generalization  of  the  singular  value 
decomposition,  and  it  solves  an  important  class  of  problems  as  well  as  exhibit¬ 
ing  their  geometric  structure. 


J.  Laurie  Snell 
Dartmouth  College 

"Random  Walks  and  Electric  Networks" 

The  connections  between  potential  theory  and  Markov  processes  are  well-known 
and  has  very  much  influenced  the  direction  of  probability  theory  in  recent 
years.  There  are  still  things  to  learn  from  these  connections.  We  show  this 
by  discussing  Peter  Doyle's  use  of  Rayleigh's  short-cut  method  to  decide  if 
discrete  random  walks  are  recurrent  or  transient.  For  this,  one  first  identi¬ 
fies  the  walk  with  an  electric  network.  Recurrence  corresponds  to  infinite 
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resistance  to  infinity  and  transience  to  finite  resistance  to  infinity.  The 
network  is  then  modified  to  obtain  a  simpler  network.  Two  kinds  of  modifica¬ 
tions  are  considered  —  shorting  and  cutting.  Sriui ting  decreases  the  effective 
resistance  to  infinity.  Thus  if  the  network,  simplified  by  shorting,  is 
recurrent  so  is  the  original.  Cutting  branches  can  only  increase  the  effective 
resistance.  Thus  if  the  network,  simplified  by  cutting  branches,  is  transient 
so  was  the  original  this  technique  is  illustrated  by  proving  recurrence  of 
simple  random  walk  in  two  dimensions  and  transience  in  three  dimensions. 


G.W.  Stewart 
University  of  Maryland 

"Matrix  Perturbation  Theory  and  Linear  Regression" 


This  talk  surveys  the  implications  of  first  order  perturbation  theory  for  the 
linear  regression  problem  with  errors  in  the  variables.  It  is  shown  how  sets 
of  regression  diagnostics  measuring  the  effects  of  these  errors  can  be  easily 
computed  from  quantities  formed  in  the  course  of  solving  regression  problems. 
It  is  also  shown  that  under  a  specific  model  for  the  errors,  the  classical 
F-tests  are  unaffected. 


B-1S 


Abstracts  of  Talks  by  Invited  Speakers  for  The  Conference  on  Applications 
of  Numerical  Analysis  and  Special  Functions  in  Statistics: 


Donald  E.  Amos 

Sandia  National  Laboratories 

"Computation  of  the  Central  and  Noncentral  F  Distributions" 


ABSTRACT 


Recursion  relations  suitable  for  rapid,  significant  digit  computation 
are  derived  for  the  cumulative  distribution  of  F'  =  (x/m)/(Y/n)  where  X  is 
X  (x,m),  Y  is  ^(n)  with  X  and  Y  independent.  The  cumulative  for  F'  is  given 
by  Bayes  theorem, 

P(F'<  f]m,n,\)  =  /°°G  (X  <  fttf/njY  =  y)p(y)dy  (l) 

0  v 

and  where  the  cumulative  noncentral  yf (\,m)  is 

Gp(x<  x)  =  |  j*  (f)P/V(x+zV2ipU/^)a z  (2) 

0  '  ' 


and  the  y8 (n)  density  is  given  by 


p(y)  * 


e"y/yi*n/2 

2n/2T(n/2) 


y20 

n  S  1 


(3) 


I  is  a  modified  Bessel  function  of  the  first  kind  with  p  =  (m-2)/2. 
If  we  integrate  (l)  by  parts  using  (2)  and  (3),  we  get 


P(F'£  <■).!-  |  (ff/2''X/2  J  to 


C1*) 
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where  a  =  fm/n  and  the  incomplete  gamma  functions  are  defined  by 
the  integrals 


/ 


-ta-l 
e  t  at 


r(a,x) 


I 


-t 

e 


t 


a-1 


dt 


a  >  0 
i^O  . 


a,  while  depending  upon  m  and  n,  is  treated  as  an  independent  variable. 


The  recursive  nature  follows  directly  by  applying  the  forward 
recursion  relation 


r(e+i,x)  _  vie,*)  *Ee~* 
rie+i)  ‘  rip)  r(p+i) 

to  (*»)  since,  with  F  the  cumulative  normal, 

r(£,x)  =  erfc(Jx)  =  zjri  F(-,JSc) 
and  T(l,x)  =  e  x 


(5) 


(6) 


,  x  &  0 

can  be  identified.  If  we  multiply  (5)  by  the  factors  of  (U)  and 
integrate  we  get 


>1  b 


e  P0  +  Q- 


(7) 


where 


(8)  . 


S  -  Mf  rT^T  f  - 


V 

and  P  «=  P(F'  <  f)  when  6  =  n/2.  In  order  to  exploit  (5),  Q-  oust 
P  P 

be  evaluated.  A  change  of  variables  in  (8) ,  z  =  t2 ,  produces  a 

form  for  Qg  which  can  be  identified  by 


(bt)dt  - - 1 - 

vx  2\h-1&u+v 


P(yF)  Lv  _bf\ 
fwoy  Vl»  Ua2/  » 


•  (9) 


a  >  0,  b  >  0  , 


ll^v  >  0,  where  $  is  the  confluent  hypergeometric  function  defined 
by  the  series 

.  (a)jj  x* 

»(a,c,x)  -  2_,  rr  <=  *  o.-v  uc 

k=0  * 


(e.)0  =  1  ,  (a)k  =  a(a+l)***(a+k-l)  =  '  , 


k  £  1 


Notice  that  if  a  is  a  negative  integer ,  the  series  reduces  to  a  poly¬ 
nomial,  while  for  c  =  a,  §(a,a,x)  =  exp(x).  Thus  with  a  =  *J(l+a)/2, 
b  =  «/a\,  =  2g+p+2  and  v  =  p  in  (9),  Qg  becomes 

«.  -  {&F  e'X£?>]  rcgr^iy  (-) 

where  p  =  Xa/[2(l+a)]  and  the  scaling  e  15  has  been  introduced  to 
eliminate  the  exponential  growth  in  i, 


$(p+B+l,p+l,p)  ~ 


Now  the  recursion  relation 


epAp+i' 

r(p+e+i) 


for  p  -»  co 


a$(a+l,c,x)  =  (2a-c+x)$(a,c,x)  +  (c-a)$(a-l,c,x)  (12) 

with  a  =  p+S+1,  c  =  p+1,  x  =  p  can  be  used  to  recur  forward  with 
(7)  and  (11),  giving 

Vl  =  Pg  +  % 

Vl  =  £V(a+1)3  *  Up+B+l)/(g+l)] 

,  (13) 

Vl  =  t(P+2B+1+o)§g  -  B5g_1]/(p+g+l) 


Vl  =  *  8+1  Ag+1 


where 


=  e  P»(lH-g+l,p+l,p)  , 

A  =  rfofi+1) 

e  \1+a/  d+a)P  r(g+i)r(p+i) 
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In  order  to  use  these  recurrence  relations,  we  must  have  P^,  $0> 
when  n  is  even  and  B^,  i^2  when  n  is  odd.  *  If  n  is  even, 

Pl*  *0*  *1  are  elementary  functions.  However,  when  n  is  odd,  is 

a  special  doubly  non-central  t-distribution  for  which  series  representa¬ 
tions  with  error  bounds  are  given.  $  ^  and  $^2  are  identified  in 
terns  of  derivatives  of  I  Bessel  functions  or  error  functions  depending 
on  whether  m  is  even  or  odd.  Computational  techniques  and  stability 
considerations  associated  with  the  recursive  computation  of  the  $  func¬ 
tions  are  also  discussed. 


A  quadrature  for  significant  digit  computation  of  P  based  on  (U) 
is  also  possible  to  cover  wide  ranges  of  parameters.  The  integrand  is 
bell  shaped  with  a  single  maximum  at  zQ.  The  idea  is  to  locate  from 
a  derivative  calculation,  estimate  the  spread  o  of  the  bell  in  terms  of 
a  fitted  normal  distribution  and  integrate  to  the  left  and  right  in  steps 
of  a  2ct  or  3d.  Since  the  integrand  is  positive,  no  losses  of  significance 
occur  due  to  subtraction.  This  procedure  takes  advantage  of  high  quality 
routines  for  the  special  functions  while  computing  only  the  dominant 
contributions  to  the  sum.  The  computation  of  zQ  is  facilitated  by  sharp 
algebraic  estimates  of  transcendental  functions  arising  from  the  derivative 
calculation. 


i 
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K.O.  Bowman" 

Union  Carbide  Corporation 

"Models  for  Approximating  the  Percentage  Points  of  Distributions" 

ABSTRACT 

Many  statistics  in  distributions  such  as  the  sample  skewness 
(/bj  =  m 3 /n 2 / 2 ) >  kurtosis  (b2  =  m4/m|),  sample  standard  deviation. 
Student's  t,  present  formidable  mathematical  problems  especially  under 
non-normality.  Even  under  normality,  only  recently  (Mulholland, 
Biometrika  1977)  has  the  null  distribution  of  /bj  for  moderate  sized 
samples  been  found,  although  acceptable  approximations  have  been  given 
(see  for  examples,  D'Agostino  and  Tietjen,  Biometrika  1973). 

Mulholland  had  followed  an  early  study  by  R.  C.  Gearly  (Biometrika 
1947),  and  used  a  careful  analysis  of  density  discontinuities  by 
reference  to  an  integral  equation  for  the  density  for  varying  sample 
size.  The  application  of  the  Gearly-Mul hoi  land  approach  even  to  a 
fairly  simple  case  such  as  the  null  distribution  of  kurtosis  appears  to 
be  completely  out  of  reach  mathematically. 

When  a  set  of  at  least  four  moments  of  a  statistic  exists,  several 
approximating  models  for  the  probability  integral  are  known. 

Outstanding  is  the  Pearson  system,  introduced  by  Karl  Pearson,  and 
extensively  studied  from  a  practical  point  of  view  under  the  leadership 
of  the  late  E.  S.  Pearson.  The  Pearson  system  (density  y)  arises  from 
the  differential  equation 

1  dy  _  -(x  +  a) 

y  dx  (Ax^  +  Bx  +  C) 


"Research  sponsored  by  the  Applied  Mathematical  Sciences  Research  Program, 
Office  of  Energy  Research,  U.S.  Department  of  Energy  under  contract 
W-7ii:)S-eng-26  with  the  Union  Carbide  Corporation. 


and  a  member  is  uniquely  determined  from  sca^e  anc*  location 

being  independently  adjusted.  In  view  of  the  many  types  which  can 
arise,  tabulation  of  percentage  points  (a  standard  set  being  0.50, 
0.75,  O.SC,  0.99,  0.999  and  corresponding  lower  levels)  for  a  limited 
domain  of  (/Bj,,B2)  has  only  been  completed  in  recent  years,  and  the 
corresponding  computer  program  makes  heavy  demand  on  small  computer 
f aci 1 ities. 

Completely  different  approaches  use  transformation  systems.  The 
Johnson  system  of  curves  (Biometrika  1949)  considers  the  mapping 
Z  =  y  +  6 f ( y )  where  Z  e  N(0,1).  Here  f(y)  =  log  y  produces  the  log¬ 
normal,  f(y)  =  log  (y/(l-y))  produces  the  S3  system,  and 
f(y)  =  sinh_1(y)  produces  the  Sy  system.  The  parameters  y,  <$  are 
determined  from  the  skewness  and  kurtosis  of  the  distribution  to  be 
approximated,  and  involve  intractable  mathematics,  especially  in  the 
case  of  S3,  which  requires  evaluation  of  four  transcendental 
integrals. 

Another  transformation  uses  T(x)  =  xSl  -  (1-x)52 
(introduced  by  Tukey,  and  used  on  empirical  data  by  Ramberg, 
Tadikamalla,  Dudewicz,  and  Mykytka,  Technometrics  1979)  where  x  is 
uniform  on  the  interval  (0,1). 

For  the  Pearson  system  we  have  introduced  approximations  for  the 
standard  percentage  points  at  some  11  levels  in  the  form  of  rational 
fractions  (iT|(b1,B2)/tt2(81,82))  of  the  degree  3  in  the  parameters 
Bj,  B2  Tor  a  domain  for  which  <  4.  A  linear  formulation  was  used 


over  sets  of  19  points  in  the  ( B j , B2 )  plane,  the  minimum  error  set 
being  chosen.  Rather  than  use  a  4th  degree  version,  we  considered  the 
segments  0  <  Bj  <  1,  and  1  <  0j  <  4,  giving  results  with  error  less 
than  0.5%. 

With  the  Johnson  system,  we  have  to  determine  as  simply  as 
possible,  the  values  of  y  and  5  from  those  of  Bi  and  02.  From  a  study 
of  contours,  it  turns  out  for  Sy  that  a  rational  fraction  in 
(B.,B2)  ^r  the  variable  {02  -  1/2(u4  +  2to2  +  3)}/b1  is  likely  to 
succeed,  oeir.g  equal  to  expCV^2).  The  analysis  of  the  S3  system 
involving  awkward  quadratures  is  more  demanding,  for  6  increases  to  <=° 
as  the  normal  point  is  approached,  whereas  y  tends  to  00  on  the  log- 
normal  line.  A  typical  model  used  is  for  example 

6  =  (32-Brl)efl(B2L-B2+log(l+31/9))f2 
where  f15f2  are  polynomial  in  ( B x , B2 )  of  degree  3.  Errors  in  the 
approximation  to  the  transformation  (which  once  found  yield  all 
percentage  points  in  terms  of  those  for  the  normal)  have  been  reduced 
to  a  few  percentages  for  the  three  domains  0  <  Bj  <  1,  1  <  01  <  4,  a:d 
4  <  Bj  <  9. 

We  have  confined  attention  to  the  problem  of  approximating 
percentage  points  for  theoretical  statistics  whose  first  few  moments 
exist;  note  that  in  many  cases  asymptotic  series  may  be  required  using 
summatory  techniques  to  establish  moments  evaluations.  We  nearly 
mention  the  corresponding  problem  for  empirical  data  which  though 
simpler  from  one  point  of  view,  since  moments  and  percentiles  are 
always  available,  is  much  more  difficult  when  inevitable  sampling 
errors  are  allowed  for. 
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James  A.  Cadzow 

Virginia  Polytechnic  Institute 

"Autoregressive-Moving  Average  Spectral  Estimation: 

A  New  Effective  Modeling  Procedure" 

ABSTRACT 

In  various  signal  processing  applications,  knowledge  of  the  spectral 
density  associated  with  a  zero  mean  weakly  stationary  random  time  series 
{ x^ }  plays  a  prominent  role.  This  spectral  density  is  formally  given  by 

Sx(u))  =  £  rx(n)e-jum  (1) 

0=-°° 

and  is  recognized  as  being  the  Fourier  transform  of  the  time  series' 
covariance  sequence 

where  E  denotes  the  expected  value  operator.  Clearly,  the  determination 
of  the  spectral  density  entails  a  complete  knowledge  of  the  infinite 
extent  covariance  sequence.  Unfortunately,  these  covariance  elements  are 
almost  never  known  in  typical  applications,  and,  one  must  therefore  resort 
to  estimation  techniques  for  determining  an  appropriately  accurate 
estimate  of  Sx(w).  This  estimate  is  generally  based  on  a  finite  set  of 
contiguous  time  series  observations  as  designated  by 


This  research  was  supported  in  part  by  the  Statistics  and  Probability 
Program  of  the  Office  of  Naval  Research  under  Contract  N0001 A-80-C-0 103 
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xl*  x2*  ’  *  *  ’  XN 

Unless  some  constraints  are  Imposed  upon  the  time  series'  fundamental 
nature,  however,  there  exists  a  basic  incompatability  in  estimating 
spectral  density  (1)  which  depends  on  the  infinite  set  of  covariance 
elements  from  the  finite  set  of  observations  (3).  This  difficulty  is 
alleviated  by  postulating  the  following  (p.q)1*’  order  ARMA  model  for  the 
time  series 


P  q 

x(n)  +  l  a.x(n-k)  =  £  b.e(n-k) 


(4) 


k=l 


k*0 


where  the  unobserved  excitation  (e(n)}  is  taken  to  be  a  zero  mean  white 
noise  series  of  variance  one.  As  proven  by  Wald,  any  continuous 
spectral  density  can  be  approximated  arbitrarily  closely  by  an  ARMA  model 
if  the  order  integers  p  and  q  arc  selected  large  enough. 

The  problem  to  be  investigated  is  then  that  of  estimating  the  , 
bj.  coefficients  of  this  ARMA  model  from  the  given  set  of  time  series 
observations  (3).  Although  there  presently  exist  procedures  for  accomplish¬ 
ing  this  task  (o.g.,  see  refs.  [1]— [5] ) ,  these  procedures  are  not  very 
effective  in  the  case  of  small  data  lengths  (i.e.,  N).  In  this  lecture, 
a  procedure  which  has  been  found  to  be  effective  for  both  short  and  long 
data  lengths  shall  be  developed.  A  brief  outline  of  the  procedure’s 
salient  features  will  now  be  given. 

The  procedure  for  estimating  the  ARMA  model  parameters  first  entails 
multiplying  both  sides  of  relationship  (4)  by  x* (n-m)  to  yield  the  "basic 
error  elements"  as  given  by 


e (m,nl  =  x(n)x*(n-m)  +  £  a  x(n-k)x*(n-m) 


(Sa? 


k=l 
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"  — - 

q 

=  £  b^e (n-k)x* (n-m)  for  m+1  n  <  N  (5b) 

k=0  q  +  1  ^  m  N  -  1  \ 

in  which  *  denotes  the  operation  of  complex  conjugation.  It  is  not  difficult 
to  show  that  if  the  time  series  is  an  ARMA  process  of  order  less  than  or 
equal  to  (p,q),  then  the  basic  error  elements  are  each  zero  mean  random 
variables  over  the  range  shown  on  the  m  and  n  variables.  With  this  in  mind, 
a  logical  selection  for  the  a^  coefficients  would  be  one  that  causes  each 
of  these  basic  error  elements  to  be  as  close  as  possible  to  their  expected 
value  of  zero.  This  objective  can  be  achieved  by  minimizing  the  following 
quadratic  functional 

f (ak)  =  e+We  (6) 

where  e  is  a  (N-m) (N-q-1)  x  1  column  vector  whose  elements  are  appropriate 
arrangements  of  the  basic  error  elements  (5),  W  is  a  nonnegative  definite 
square  matrix^and,  t  denotes  the  operation  of  complex  conjugate  trans¬ 
position.  This  criterion  is  seen  to  be  a  quadratic  function  of  the 
coefficients  through  the  basic  error  element  relationship  (5a).  Once  the 
optimal  set  of  a^  coefficients  have  been  thus  determined,  the  modified 
Welch  method  [6]  may  be  applied  to  identify  the  bk  coefficients  effects 
on  the  spectral  estimate. 

In  this  lecture,  a  more  detailed  development  of  this  new  ARMA  model 
method  shall  be  given.  Moreover,  the  new  method's  performance  will  be 
empirically  compared  to  such  classical  spectral  estimation  techniques  as 
(i)  the  Box-Jenkins  ARMA  method,  (ii)  Burg’s  maximum  entropy  AR  method, 
and,  (iii)  the  Periodogram.  In  this  comparison,  it  is  found  that  the  new 
method  significantly  outperforms  the  classical  methods. 


If  this  new  method  is  to  achieve  its  full  potential,  however,  a  number  of 
numerically  based  issues  need  further  investigation.  In  keeping  with  the 
spirit  of  this  conference,  these  issues  will  be  dwelled  upon  and 
suggestions  solicited.  Perhaps  the  most  significant  issue  that  needs 
further  attention  is  that  of  selecting  the  weighting  matrix  W  in 
criterion  (6).  Preliminary  empirical  evidence  attests  to  the  significance 
of  this  choice  (7}. 
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George  Cybenko 
Tufts  University 

"The  Efficient  Solution  by  Orthogonalization  of  Linear 
Prediction  Problems  for  Stationary  Time  Series" 


Abstract 

The  linear  prediction  problem  for  stationary  time  series  has 

traditionally  been  solved  by  forming  the  normal  equations  and 

solving  them  by  either  classical  or  fast  Toeplitz  algorithms. 

The  main  obstacle  to  using  orthogonalization  has  been  that  that 

2 

approach  requires  an  order  of  magnitude  more  computations  (0(Np  ) 
as  opposed  to  O(Np)). 

2 

An  0 (Np  )  orthogonalization  technique  is  described  for 
general  matrix  orthogonalization  which  yields  an  O(Np)  method  in 
the  special  case  of  linear  prediction.  Orthogonalization  is 
thereby  made  competative  with  the  normal  equations  approach. 
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The  linear  prediction  problem  for  a  stationary  time  series 
may  be  formulated  as  a  matrix  linear  least-squares  problem,  where 
the  design  matrix  has  a  Toeplitz  structure. 

Specifically,  let  x(n)  be  a  time  series  with  only 
finitely  many  nonzero  terms,  so  that  x(n)=0  for  0<n<N  say.  The 
linear  prediction  problem  of  order  p  is  to  find  coefficients 
a.,..., a  minimizing  the  expression 

**■  r 


■  Z 


(  x (n) 


+  a^  x(n-l)  +  . 


+  ap  x(n-p) ) 


where  the  summation  is  over  all  n. 


This  problem  occurs  in  a  variety  of  applications:  Wiener 
filtering,  stochastic  model  identification,  speech  analysis  and 
synthesis,  and  geophysical  signal  processing,  to  name  a  few. 


Letting 


a  = 


x(afl 


x  (N) 
0 


0 
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0 

0 

0 

x(0) 

0 

0 

x(l) 

x{0) 

0 

x  (2) 

x(l) 

x  (0) 

0 

0 

0 

0 


0  0  0 
0  0 
0  0 


.  x (N-2) 
X (N)  x(N-l) 
0  x  (N) 


the  linear  prediction  problem  is  equivalent  to  solving  the  matrix 
equation 


V  ■  6 


in  the  least-squares  sense. 

The  traditional  approach  to  solving  these  equations  has  been 
to  form  the  normal  or  Yule-Walker  equations 
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T 

XpXp  a 


x 


which  is  a  simple  system  of  equations.  Because  XpXp  has  a 

Toeplitz  structure  also,  it  can  be  computed  in  O(Np)  operations 

3  2 

after  which  the  pxp  system  requires  0(p  )  or  0(p  )  more  for 

solution.  Since  N  is  always  much  larger  than  p,  the  O(Np)  term 

dominates.  Now  it  is  well  known  that  for  small  residuals, 

solving  least-squares  problems  by  orthogonization  is  more 

accurate  than  by  forming  and  solving  the  normal  equations. 

However,  the  orthogonalization  of  Xp  by  any  of  the  classical 

2 

methods  requires  0(Np  )  ,  or  an  order  of  p  more,  steps  than  the 
previously  described  approach.  For  this  reason, 

orthogonalization  has  not  been  used  to  solve  linear  prediction 
problems. 

In  this  paper,  we  present  an  orthogonalization  technique 

2 

which  appears  to  be  new,  and  which  is  0(Np  )  for  general 
unstructured  matrices,  and  streamlines  to  an  O(Np)  method  for  the 
special  Toeplitz  structures  that  arise  in  linear  Prediction. 

It  is  shown  that  this  procedure  is  the  same  as  the 
Itakura-Saito-Burg  "lattice  method"  for  linear  predictive 
filtering  and  deconvolution. 

Although  a  complete  error  analysis  is  not  presently 
available,  partial  results  indicate  that  the  algorithm  has  good 
accuracy  properties. 
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Robert  B.  Davies 

University  of  California,  Berkeley 
"Maximum  Likelihood  Estimation" 


ABSTRACT 

Estimation  of  parameters  for  "non-standard"  statistical 
problems  frequently  involves  finding  the  maximum  of  the 
likelihood  function  or  some  other  function  related  to  it. 
Usually  analytic  solutions  won^t  be  possible  and  a  numerical 
method  will  be  required.  While  a  variety  of  optimizing  pro¬ 
grams,  available  in  the  standard  subroutine  libraries  can  be* 
used,  the  likelihood  function  has  particular  properties 
which  it  is  worth  taking  advantage  of  and  at  the  same  time 
the  statistician  has  particular  requirements  which  the  stan¬ 
dard  optimization  programs  do  not  necessarily  provide. 

This  paper  concentrates  on  the  situation  where  the  likeli¬ 
hood  function,  its  first  derivative  and  an  approximation  to 
the  second  derivative  are  available.  Suppose  8  denotes  the 
vector  of  unknown  parameters;  L(0)  the  log-likelihood;  A(0) 
=  dL(0)/d(0),  T(0)  =  EQ[A(0)  A'(0) ]  =  -EQ[d2L(0)/d02] .  I  will 
suppose  -T(0)  is  used  as  the  approximation  to  the  matrix  of 
second  derivatives. 

In  some  problems,  once  L  has  been  calculated,  not  too  much 
extra  work  is  required  to  obtain  A  and  Tor  if  L  and  A  are 
calculated  very  little  work  is  required  to  find  r  as  well. 
In  these  situations,  it  seems  especially  worthwhile  using 


optimization  methods  that  involve  L,  A,  and  r  if  the  calcu¬ 
lation  of  L  is  expensive,  if  the  number  of  parameters  is 
large,  or  if  the  problem  is  ill-conditioned.  Note  that  P(9) 
is  non-negative  definite  (and  its  being  singular  usually 
indicates  an  error)  so  that  one's  optimize  program  need  not 
take  into  account  the  possibility  of  r(8)  having  negative  or 
zero  eigenvalues.  The  basic  iterative  step  is 

Vi  -  ®k  +  r'1(ek>a(ek> 

but  we  have  found  that  some  kind  of  line  search,  based  on  L 
and  A,  is  advisable. 

According  to  maximum-likelihood  theory  the  variance  / 
covariance  matrix  of  one's  estimates  is  approximately  r  ^(9) 
so  one  would  normally  want  this  printed  out,  or  rather,  one 
wants  the  standard  errors  and  correlations  derived  from  it. 

Once  one  has  fitted  a  model  one  will  frequently  want  to  test 
whether  it  should  be  extended  by  fitting  further  parameters. 
One  possibility  is  simply  to  program  the  extended  model,  fit 
it,  look  at  the  increase  in  the  log-likelihood  and  so  carry 
out  a  likelihood  ratio  test.  But  a  simpler  approach  is  to 
carry  out  Neyman's  C(<x)  test.  For  this  one  still  needs  to 
calculate  A  and  P  for  the  full  model  but  only  for  values  of 
9  corresponding  to  the  old  model.  One  then  calculates  A' r  *A 
at  the  values  of  9  corresponding  to  the  old  model.  If  the 
extended  model  is,  in  fact,  not  necessary  this  term  will 


have  an  approximately  chi-squared  distribution. 
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This  also  provides  a  very  convenient  stopping  criterion  for 
the  iterative  fitting  process.  We  test  the  hypothesis  that 
our  current  value  of  0  is  correct  by  calculating  A'  f1  ^A  at 
each  step.  However  rather  than  using  the  critical  points  of 
the  chi-squared  distribution  for  deciding  when  to  stop,  one 
stops  when  A'  P  <  .001  say.  Another  way  of  looking  at  this 
would  be  to  say  that  A'r~^A/2  approximates  the  possible 
improvement  in  the  likelihood.  In  any  case  the  objective  is 
to  make  sure  the  difference  between  each  estimate  obtained 
and  the  actual  maximum-likelihood  estimate  is  small  compared 
with  the  standard  error  of  the  estimate. 

If  one  does  calculate  A'r-1A  at  each  step  then  the  program 
can  also  be  used  for  carrying  out  a  C(«)  test  if  it  can  be 
run  for  just  one  step. 


Sometimes  A  and  r  naturally  factorize:  A  =  X'Y,  P  =  X'X 
where  now,  perhaps,  the  components  of  Y  are  closely  related 
to  the  original  observations.  In  this  case  one  can  use  a 

t 

Householder  transformation 


X  =  Q 


R 
0  /' 


Y  =  Q 


Q  orthogonal  so 

that  the 

9k+l  " 

ek  +  r-1zi; 

the  stop- 

.001} 

and  P  1  = 

(R'R)"1. 

This  approach  can  be  important  when  T  is  ill-conditioned 
and,  in  fact,  there  seems  to  be  more  justification  for  it  in 
the  present  context  than  in  the  usual  linear  least  squares 


i 


situation. 


This  Householder  approach  is  also  particularly  convenient 
when  one  wants  to  use  recursive  fitting  techniques.  That  is 
when  9  is  broken  into  two  (or  more)  parts; 


and  for  each  iteration  of  ©2  one  selects  9^  to  maximize  the 
likelihood. 

Considering  the  selection  of  the  initial  values  from  which 
to  start  the  interative  process.  The  techniques  we  have 
considered  so  far  may  not  be  effective  if  one  is  a  long  way 
from  the  correct  value.  Global  searching  may  be  appropriate. 
However  we  may  be  able  to  begin  model  fitting  with  a  very 
simple  model  involving  only  one  or  two  parameters  which  may 
be  able  to  be  estimated  by  other  means.  In  particular,  all 
"treatment  effect"  parameters  can  be  set  to  zero.  After  fit¬ 
ting  the  simple  model  other  parameters  can  then  be  added  and 
fitted  until  the  full  model  is  fitted.  One  important  attri¬ 
bute  of  the  optimize  program  would  be  the  ability  to  fix 
certain  parameters  and  temporarily  eliminate  them  from  the 
fitting  process. 

Our  comments  have  been  specifically  for  the  likelihood  func¬ 
tion.  In  some  statistical  problems  it  is  appropriate  to  max¬ 
imize  some  other  function,  W<0)  say.  A  maximum  likelihood 
program  can  be  applied  identically  to  the  maximization  of 

such  functions  which  satisfy  the  usual  regularity  conditions 
and  also  EQ[dW/d9]  *  0,  covQ [dw/d9,d (L-W) /d9)  =  0. 
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Charles  F.  Dunkl 
University  of  Virginia 

"Discrete  Orthogonal  Polynomials" 

ABSTRACT 

The  binomial  and  hypergeometric  are  important  discrete  distributions. 
The  orthogonal  polynomials  for  the  one-variable  case  are  the  well-known 
Krawtchouk  and  Hahn  polynomials .  By  considering  some  several -variable 
analogues  one  is  led  to  another  family  of  hypergeometric  polynomials 
expressed  as  balanced  ^F^-series  (appearing  in  physics  as  6-j  symbols) . 
Special  cases  of  these  give  new  representations  of  Krawtchouk  and  Hahn 
polynomials  associated  to  symmetric  distributions  (that  is,  binomial  with 
p  =  h,  hypergeometric  * 


1.  The  Hahn  polynomials: 

For  integer  parameters  a,b,c  we  have  the  hypergeometric  distribution 
0  (c-x)  ’  max(0>c-k)  x  £min(a,c).  The  family  of  orthogonal  poly¬ 

nomials  for  this  weight  is  E  (a,b,c,x)  =  T™  (-l)"1  (?)  (b-crt-1)  .  (a-m+1) 

m  Lj=0  j  j  m-j 


(-x).(x-c)  .  =  (-l)m(-a)  (-c)  ,F  (  ra’m  3  b  X  ;1) ,  0  <  m  <  min(a,b,c,a+b-c) 

j  m-j  ra  m  j  L  -a ,  -c  _  — 


This  is  a  useful  notation  for  the  Hahn  polynomials  (see  Karlin  and  McGregor, 


Scripta  Math.  26  (1961),  33-46)  and  the  relation  is 

Em(a,b,c,x)  =  (-l)m(-a)m(-c)raQm(x;-a-l,-b-l,c) 

(the  present  notation  exhibits  certain  symmetries,  and  is  polynomial  in  the 
parameters).  These  functions  appear  as  intertwining  functions  on  the  symmetric 
group  (consider  the  set  of  2x2  contingency  tables  with  row  totals  a  and  b,  and 
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column  totals  c  and  a+b-c;  for  any  permissible  arrangement  of  a+b+c  objects 
in  a  2x2  array,  there  is  a  subgroup  of  permutations  preserving  rows,  and 
another  preserving  columns).  This  fact  has  been  used  by  C.D.  (SIAM  J.  Math. 
Anal.  9  (1978),  627-637)  to  derive  product  and  addition  (b  =  c  only)  formulas. 


2.  Hahn  polynomials  in  several  variables: 


For  positive  integers  r,  . . .  ,ln  we  have  the  hypergeometric 


ft  a  12  n 

distribution!  ]...(  \  with  T.x.  =  r. 

VX1 J  \\  ) 


is  given  by 


A  set  of  orthogonal  polynomials 


n-1  j-1  n-1 

n.  .  E  (£.,Jn  ...  SL.-2  m  ,  r  -  £  x  -  £  m.,x.) 

j=l  j  ^i-j+1  i  \=j+1i  i  i£j+1  1  J 


where  m^.m^, . . . are  nonnegative  integers  subject  to  various  constraints 

(see  §1).  The  product  is  a  polynomial  in  x, ,...,x  ..of  degree  1. m. 

i  n-l  i  i 

2 

(Pearson's  X  -statistic  for  the  2*n  contingency  table  can  be  expressed  as  a 
weighted  sum  of  squares  of  the  linear  polynomials) .  These  polynomials  were 
studied  by  Karlin  and  McGregor  (pp.  261-288  in  "Theory  and  Application  of 
Special  Functions",  R.  Askey,  ed.,  Academic  Press  1975)  in  the  Q-form  for 
a  genetics  application,  with  the  i^'s  being  replaced  by  negative  real  numbers. 
The  present  situation  is  associated  to  the  permutation  groups  of  a  2*n 
contingency  table  (column  totals  row  totals  r  and  £  i^-r)  .  The 

given  basis  depends  on  the  ordering  of  the  variables,  thus  a  rearrangement 
will  produce  a  different  basis.  The  connection  coefficients  between  two  of 
these  bases  form  a  set  of  orthogonal  transformations,  one  for  each  total  degree 
I^m^.  In  the  case  n  =  3,  given  the  degree  N,  there  is  one  free  parameter  m^ 
(note  m^  =  N-m^) ,  so  the  connection  coefficients  form  a  family  of  orthogonal 
functions  in  one  variable.  These  were  determined  to  be  (C.D.,  Pacific  J.  Math 
to  appear)  balanced  ^F ^-polynomials,  and  will  be  discussed  in  the  next  section. 


B-39 


3.  The  discrete  ^F^-polynomials: 

For  an  integer  N  and  parameters  a,b,c  (positive  integers  in  the 
contingency  table  or  group  case),  define  a  weight  function 


W(x 


K  ,  ( b+c-2x+l\  (-a)N_x(-Ox(-N)x(N-a-b-c-l)x(-l) 

,a’  ,C  \b+c-x+l  /  x! (-b)x(x-b-c>N 


(where  x  =  0,1,..., N  and  a-N  £  x  <_  min(b,c,b+c-N)  if  a,b,c  are  integers) 
and  the  hypergeometric  polynomial 

/-k,k-a-c-l,-x,x-b-c-l  A 
pk(x;a,b,c)  =  4F3  ^  -N , -c , -a-b-c+N-1  ;1J 

(a  balanced  series;  sum  of  denominator  parameters  exceeds  sum  of  numerator 

2 

parameters  by  1) ,  a  polynomial  of  degree  k  in  (x-(b+c+l)/2)  .  The 
orthogonality  relation  is 

£xw(x;a,b,c)pk(x;a,b,c)p£(x;a,b,c)  =  6RJl/w(k;b,a,c)  . 

J.  Wilson  (SIAM  J.  Math.  Anal.  11  (1980),  690-701)  showed  that  the  Racah 
6-j  symbols  could  be  expressed  as  ^F^-polynomials  (his  parameters  are  different, 
and  he  also  found  continuous  orthogonality  relations,  N  replaced  by  a  real 
number) . 


4.  Special  cases  of  the  ^F^-polynomials: 

Recognizable  and  interesting  distributions  can  be  obtained  by  taking 

b+c+1  to  be  an  integer,  say  -s,  with  s  ^  0.  Neglecting  constants  one  obtains 

(-c)x 

the  weight  (N+s+xHN!x)(c+sii)  A(S’X)’  X  =  0,1, ...  ,N  ^ere 

x 

A(0,x)  =  (2-6  ),  and  for  s  >_  1, 

XU 

A(s,x)  =  (x+l)g  ^(2x+s)/s!,  the  weight  being  positive  for  -s-1  <  c  <  0. 
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2 

In  general  one  obtains  polynomials  not  in  x,  but  rather  in  (x+s/2)  .  However 

in  special  cases  (s  =  0,  -1  <  c  <  0  or  s  >  1,  c  =  -s/2)  one  can  piece 

together  two  ^F^-families  to  get  one  family  of  polynomials  in  x,  a  situation 

analogous  to  the  quadratic  transformation  for  Legendre  polynomials,  relating 

P„  (x)  to  P^’  ^ (2x2-l)  and  P„  ,,(x)  to  xP^’^  (2x^-1) .  The  symmetric  Hahn 
Zn  n  2n+l  n 

3  3 

polynomials  (weight  (^) (  ^))  form  an  example  of  this.  Here  is  the  even  case 

(c  =  2N)  (the  odd  case  c  =  2N+1  involves  s  =  1)  : 


Q2n(x;-a-l,-a-l,2N)  = 
Q2n+1(x;-a-1’-a-1’2N) 


(%)n(N-a)n 

H^n 


4F3 


(-n,n-a-4,N-x,x-N  \ 

-N,5s,N-a  ;1J 


/N-x\  (3/2)n(N~a+1) 

V  N  )  (5s-N)  (-a) 

x  #  n  n 


(-n , n-a+^.N-x+l , x-N+1 
-N+1,3/2 ,N-a+l 


Let  a  =  -1  to  get  the  discrete  Chebyshev  polynomials  (discrete  uniform 

4N 

distribution)  ,  a  =  2N-%  for  the  weight  (2x) ,  or  let  a  to  get  the 

symmetric  (p  =  h)  Krawtchouk  polynomials,  for  example 

(h)  (~ n,N-x,x-N  \ 

K2n(x;li,2I0  =  Tia-K) n  3F2  ^  -N,%  ;1J  " 


- 
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Roy  S.  Freedman 
Hazeltine  Corporation 

"Random  Walks  and  Statistical  Communication" 


ABSTRACT 


1.  INTRODUCTION 


We  discuss  an  application  of  special  function  techniques  to 
a  problem  that  occurs  in  the  analysis  of  bit  error  probabilities  in 
a  certain  digital  communication  system.  Let 


(  U/  - 


«  -  1_ 


) 


^Cosg,- 

i*  I 


-  /  A  ( 

k  -  I  51 


fj  (  H  )  —  Pi~o  £>  (  L  > 

Pi  (**i>*i)  —  P‘0  4  ( -Si 


In  the  above  expressions,  the  j§;  are  independent  and  uniformly 
distributed  random  variables  on  [  ;  n,  m,  w  are  positive 

integers  with  n,  m  w  and  j  =  /\f~~ l"«  The  problem  is  to  determine 

P!  (m'  n>  and  ?2  (m,  n) . 


These  probabilities  can  be  interpreted  in  the  following  sense. 

We  are  given  two  drunks:  each  drunk  takes  a  unit  step  in  an  arbitrary 
direction.  The  origin  is  fixed  at  the  point  where  drunk  B  starts; 
drunk  A  is  given  a  head  start  of  (w-n)  steps  in  the  x  -  direction. 

After  a  certain  time,  drunk  B  takes  n  random  steps  and  drunk  A  takes 
m  random  steps  (so  the  total  number  of  steps  that  drunk  A  takes  is  w) . 

(m,  n)  is  the  probability  that  drunk  B  is  farther  from  the  origin 
than  drunk  A.  In  this  case,  "farther  from"  is  interpreted  to  mean 
that  the  resultant  magnitude  of  drunk  B  is  greater  than  the  resultant 
magnitude  of  drunk  A.  If  we  are  only  interested  in  how  far  the  drunks 
have  traveled  in  a  certain  direction  (without  loss  of  generality,  the 
x  -  direction} ,  then  the  probability  that  the  sum  of  the  x  -  components 
of  drunk  B  is  greater  than  the  sum  of  the  x  -  components  of  a  drunk  A 
is  given  by  P2  (m,  n) .  A  discussion  of  the  original  interpretation 
in  the  communication  system  investigated  is  given  in  Section  4.  It 
suffices  at  this  point  to  denote  P^  (m,  n]  as  the  probability  for  the 
non-coherent  case  and  P2  (m,  n)  as  the  probability  for  the  coherent 
case.  Expressions  for  P^  (m,  n)  and  P2  (m,  n)  will  be  derived  below. 
Sinplificatians  involve  the  use  of  Bessel  function  identities  as  well 
as  the  use  of  generalized  functions. 


It  is  well  known  that  the  characteristic 


2 .  NON-COHERENT  CASE 
Let  . 

function  of  X  in  polar  coordinates  (see  reference  CU)w 

■) 

The  addition  of  independent  random  variables  yields  (with  k  =  w-m) 
&(*)  = 

The  first  factor  is  the  characteristic  function  of  the  sum  of  m 
random  phasors.  The  second  factor  is  the  characteristic  function 
of  the  "head-start."  We  note  that  the  head-start  can  be  taken  in  an 
arbitrary  direction  of  length  k  if  all  the  other  random  angles  are 
taken  with  respect  to  the  first  head-start  angle  as  a  reference. 

A  similar  argument  shows 

<Ps(t)  =  7C * 

The  density  and  distribution  functions  of  a  and  b  follow  from  Hankel 
Transformations  and  suitable  Bessel  function  identities.  Consequently, 
an  expression  for  Cm,  n]  is 

o  o  O  . 

The  innermost  integral  (with  respect  to  u)  is  the  density  function 

of  b.  The  middle  integral  (with  respect  to  t)  is  the  cumulative 
distribution  function  of  a;  i.e. ,  the  probability  that  a  is  less 
than  y. 

After  rearranging  terms  we  note  that 

fV’t  7.1]*)^  =■  Si*-*) 

*0 


This  formula  can  be  proved  by  using  the  discontinuous  integral  of 
Weber  and  Schafheitlein.  In  reference  (2)  it  is  equation  (29)  on 
page  51.  There  seems  to  be  a  typographical  error  in  the  convergence 
criteria  in  this  equation.  It  should  read 

$e.  ['O+sti.'-f  +  l)  y° 

The  parameters  of  interest  are  _f  ~ ~  0  and  7?  •“  ^ 

In  this  case,  the  hypergeometric  function  is  completely  degenerate 
and  the  Weber-Schafheitlein  integral  reduces  to  the  step  function 
represented  in  reference  (3)  as  equation  (9),  page  406.  If  we 
differentiate  this  expression  twice,  we  obtain  the  unit  doublet  as 
indicated  above.  The  integral  for  (m,  n)  simplifies  to 

This  use  of  the  Weber-Schafheitlein  integral  seems  to  have  first 
been  used  by  Doyle  (see  reference  (4)).  Kluyver  uses  the  integral 
to  derive  an  expression  for  the  density  function  of  b  (see  ref¬ 
erence  (3),  page  419). 

After  an  integration  by  parts  we  obtain 

f  r  &  ^  1 

(*,«)  =  Ji_  1"  Kj  J  («t)  J , 'Vo  A* 


The  last  integral  cannot  be  simplified  further.  For  n  and  m  large 
we  can  obtain  asymptotic  expressions  for  the  density  of  a  and  b: 
b  is  asymptotically  Rayleigh  distributed  and  a  is  asymptotically 


Rice-Nakagami  distributed.  We  find  that 

—  ^  /(  m  ) 


PjK*)  = 


*1 


d 


a  result  which  agrees  exactly  for  k=0. 

3 .  COHERENT  CASE 

In  this  case,  we  do  not  deal  with  magnitudes  of  phasors  and 
do  not  calculate  densities  with  Hankel  Transforms.  Consequently, 
the  following  characteristic  functions  for  A  and  B  are 

X  X 

A* 

<P.Lt)= 

After  taking  Fourier  Transforms  and  making  use  of  the  identity 

-L  f  Cos 

7 r  .1 


we  find  that 


♦  M  ,  & 


P L  Lm*  *  )  ~  J'JT' 


U 


-  CO  <J 


CoS 


t(K-K)  X 


f*  r & 


II 


— oO  °  Vo 


cs  t  (*.-<)  *1  ix  it 

y 


which  after  simplification  yields 
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fi  C*,*) 


X  -  -L.  f  J*. 

oC  ^  \  0  %- 


For  n  and  m  large,  a  and  b  are  respectively  asymptotically  non-zero 
and  zero  mean  gaussian  distributed.  Consequently 


u,tk  crx-  <b\  +  M/sl  • 


This  result  also  agrees  with  the  exact  expression  for  k=0 , 

4.  INTERPRETATION  FOR  STATISTICAL  COMMUNICATION 

In  our  system  either  a  "mar k"  or  "space"  signal  is  transmitted, 
both  events  being  equally  likely.  We  assume  a  mark  is  transmitted. 
The  probability  of  error  is  given  by  Prob  (receive  space  given  that 
a  mark  was  transmitted) .  This  probability  depends  not  only  on  the 
transmission  media  but  also  on  the  encoding  and  modulation  schemes . 
We  will  assume  that  a  mark  (or  space)  is  encoded  onto  a  set  of  w 
waveforms,  with  each  waveform  suitably  represented  as  a  phasor.  We 
will  assume  that  a  phase  reference  is  known  at  the  receiver.  The 
effect  of  noise  is  to  change  the  amplitude  and  phase  of  a  particular 
waveform.  Intuitively,  the  more  waveforms,  the  less  likely  it  will 
be  that  noise  will  disrupt  all  w  waveforms.  At  the  receiver,  each 
particular  waveform  will  be  detected  and  "hard-limited"  in  that  the 
magnitude  of  a  received  phasor  will  be  normalized  to  unity  if  that 
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magnitude  exceeds  a  certain  threshold.  This  minimizes  the  effects 
of  strong  interference;  information  is  still  contained  in  the  phase. 
We  assume  that  in  the  w  mark  waveforms,  m  contain  interference  and 
in  the  w  space  waveforms,  n  contain  interference.  In  the  non¬ 
coherent  case,  a  mark  is  decided  if  the  magnitude  of  the  resultant 
phasor  sum  for  the  mark  waveforms  is  greater  than  the  resultant 
phasor  sum  for  the  space  waveforms.  In  the  coherent  case,  we 
compare  the  sum  of  the  values  of  the  components  in  each  waveform. 

The  error  probabilities  are  then  given  by  (m,  n)  and  ?2  (m,  n) . 
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"A  Stochastic-Dynamic  Model  for  Global  Atmospheric 
Mass  Field  Statistics" 


ABSTRACT 


A  model  that  yields  the  spatial  correlation  of  atmospheric  temperature  data 
has  been  developed.  It  involves  the  solution  of  the  potential  vorticity  equation 
forced  by  random  noise:  ' 

(v C„  e)  $ (\ , V  ( 0;  u,)  (1 ) 


2 

where  V  is  the  laplacian  operator  in  the  unit  sphere,  X  and  ©are  longitude 
and  latitude,  <^>  is  the  temperature  and  F  is  white  noise  corresponding  to  a 
random  realization  6l)  • 

The  spatial  correlation  P  is  then  computed  from 


f(X  e,  3  v , et)  =  ‘2> 


where  E  is  the  expected  value. 


Three  methods  of  solution  have  been  tested.  In  the  first  method,  Eq.  (1)  was 
solved  by  expansion  in  spherical  harmonics  and  the  correlation  function  was 
computed  analytically  using  the  expansion  coefficients.  In  the  second  method, 
the  finite  difference  equivalent  of  Eq.  (1)  was  solved  using  a  Fast  Poisson 
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Solver.  The  correlation  function  I  was  computed  using  stratified  sampling  of 
the  individual  realizations  of  F(zo)  and  hence  of  ^(co).  In  the  third  method 
an  equation  for  F  was  derived  from  Eq.  (1)  and  solved  directly  in  finite 
differences  by  succesive  applications  of  the  Fast  Poisson  Solver.  The  three 
methods  were  compared  for  accuracy  and  efficiency,  and  the  third  method  was 
chosen  as  clearly  superior. 

The  results  agree  well  with  the  latitude  dependance  of  observed  atmospheric 
correlation  data.  The  value  of  the  parameter  C0,  chosen  by  best  fit  to  the 
data,  is  close  to  the  value  expected  from  dynamical  considerations. 

These  results  provide  the  basis  for  an  optimal  choice  of  coefficients  for 
statistical  analysis  of  atmospheric  data.  , 
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ABSTRACT 


1.  Summary.  The  concept  of  Standardized  Generalized  Variances  (SGVs) 
is  introduced.  Several  new  problems  of  multivariate  statistical 
inference  are  formulated  on  the  basis  of  these  SGVs.  It  is  shown 
that  in  addition  to  providing  several  new  statistical  tests,  many 
existing  problems  of  multivariate  tests  of  significance  can  be 
regarded  as  special  cases  of  these  formulations  and  can  also  be 
extended  to  their  full  generalities.  The  null  and  non-null  distribu¬ 
tions  of  these  test  criteria  are  thus  of  vital  importance.  Consider¬ 
ing  multivariate  normal  populations  with  general  and  equi-correlated 
covariance  matrices,  these  distributions  are  deduced  in  computable 
forms  in  terms  of  a  variety  of  Special  Functions,  e.g.,  Pincherle's 
H-function,  Meijer's  G-function,  Rummer's  function,  Whittaker's 
function,  Riemann's  Zeta  function  and  Psi-function.  The  highly 
desirable  property  of  unbiasedness  is  also  established  for  most  of  the 
above  test  criteria.  Finally,  applications  of  the  above  tests  to  a 
wide  spectrum  of  applied  research  are  also  illustrated  by  examples 
taken  from  the  existing  literature. 


B-51 


2.  Introduction.  Let  X  be  a  p-dimensional  random  vector  variable 
with  dispersion  matrix  Z.  Two  well-known  measures  of  multidimensional 
scatter,  obtained  by  generalizing  the  variance,  the  univariate 
measure,  are  Z  and  the  generalized  variance  (GV) ,  |z[  =  det(E),  intro¬ 
duced  by  Wilks  (1932).  For  multivariate  normal  populations,  Likeli¬ 
hood  Ratio  Tests  (LRTs)  for  £s,  of  course  of  same  dimensionalities, 
and  some  optimum  properties  of  these  tests  are  known.  But,  when 
multi-dimensional  scatter  of  populations  of  different  dimensions  need 
to  be  compared,  these  tests  cannot  be  defined.  However,  using 
which  we  will  nomenclature  as  Standardized  Generalized  Variance  (SGV) , 
such  comparisons  become  meaningful.  Since  |z|  represents  the  volume 
in  p-d intensions,  note  that  |l|^p  becomes  a  measure  so  scaled  as  to 
become  comparable  with  scatter  for  a  scaler  random  variable.  Apart 
from  this  generality,  need  for  tests  of  generalized  variances  have 
been  also  felt,  on  its  own  right.  [i|,  being  a  scalar,  is  more  suit¬ 
able  and  easier  to  work  with  than  the  matrix  Z.  Hoel  (1937)  was 
probably  first  to  realize  this  need  and  later  Eaton  (1967)  studied 
some  problems  of  statistical  inference  associated  with  a  single  GV. 

The  GV  has  been  extensively  used  in  applied  research,  e.g.,  by 
Goodman  (1966)  in  Agricultural  Statistics,  Gnanadesikan  and  Gupta 
(1970)  in  Ranking  and  Selection,  Arvanitis  and  Afonja  (1971)  in  Sample 
Survey,  etc.  While  the  estimation,  e.g.,  van  der  Vaart  (1965), 
Shorrock  and  Zidek  (1976),  and  the  distribution,  e.g.,  Bagai  (1965), 
Mathai  (1972)  have  been  studied  in  some  detail,  little  seems  to  be 
known  about  tests  for  GVs.  This  paper  attempts  to  bridge  that  gap 
via  the  extensive  use  of  Special  Functions. 
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Suppose  Xj,s  are  independently  distributed  as  N 


being  general  dispersion  matrices,  i=l,...,k.  LRTs  are  deriv  i  for 


1/p.  2  '-/p. 

Hqi  :  =  CJq  (given)  >  0,  for  some  fixed  i;  HQ2  :  JzJ  1  = 


l/P* 


Z. |  ^ ,  for  some  fixed  i  and  j  and  finally  H  ,  :  jl, 


.  !/P,- 


J'  -  '  03  •  *“x 

i=l,...,k  against  appropriate  two-sided  alternatives.  The  test 


all  equal. 


criteria  turn  out  to  be  quite  elegant  multivariate  analogues  to  those 

in  the  univariate  cases.  The  special  case,  when  s  are  equi- 

correlated  matrices,  i.e.,  Z.  «*  (.p  )  =  £  ,  where  .p  =  1,  u  =  v. 

i  x  uv  p .  x  uv  ’  ’ 


and  . p  =  p.,  u  ^  v,  is  next  considered.  As  Anderson  (1963)  has 
x  uv  1 


pointed  out,  the  statistical  inference  dealing  with  correlation 
matrices  can  become  much  more  complicated  compared  to  those  dealing 
with  covariance  matrices.  However,  such  equi-correlated  structure  as 
considered  above  is  of  extreme  importance  [see,  e.g.,  Kshirsagar 
(1978),  p.  227]  and  has  extensive  applications  [see,  e.g.,  Mitra  and 
Ling  (1979)]  in  applied  research.  After  exhibiting  the  shortcomings 
of  the  LRT  in  this  case,  some  new  tests,  including  one  based  on  the 
smallest  chracteristic  root  of  the  dispersion  matrix,  are  given  in 
this  paper. 

The  solutions  to  the  distributional  problems  associated  with 
the  various  test  statistics  considered  above  need  extensive  use  of 
Special  Functions.  The  exact  distributions  for  both  the  null  and 
non-null  cases  are  presented  for  most  of  the  above  test  criteria. 

The  percentage  points  of  these  distributions  can  be  obtained  from 
existing  mathematical  tables  since  the  distributions  are  represented 
in  suitable  computable  forms.  Examples  of  construction  of  tables  and 
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the  general  procedure  of  obtaining  them  for  such  computable  forms  of 
the  distributions  exist  in  current  literature,  e.g.,  Mathai  (1979). 
Further,  many  of  the  existing  tables  can  also  be  exploited  to  give 
the  percentage  points.  Large  sample  approximations  to  the  above 
distributions  are  also  presented. 

Considering  general  of  equal  dimensions,  for  and 
it  is  shown  that  the  'modified'  LRTs  are  unbiased — a  result  parallel 
to  Sugiura  and  Nagao  (1968)  on  tests  of  covariance  matrices.  Total 
and  partial  unbiasedness  of  some  of  the  tests  proposed  for  the  equi- 
correlated  case  are  also  established.  It  is  shown  that  the  same 
invariant  measure  under  the  full  linear  group  that  was  exploited  to 
prove  the  unbiasedness  in  the  cases  of  tests  for  covariance  matrices 
can  also  be  exploited  here  to  establish  the  unbiasedness  of  tests  for 
SGVs. 


3.  Applications.  In  addition  to  the  mathematically  interesting 
nature  of  the  problem  and  the  applications  cited  above,  there  lies  a 
rich  fertile  area  for  numerous  applications  of  the  SGVs.  In  fact, 
wherever  variance  is  employed  for  univariate  situations,  SGVs  seem  to 
be  applicable  for  the  multivariate  situations.  Some  examples  are 
cited  below. 

(a)  Multivariate  Quality  Control.  It  is  well  known  [e.g.,  see 
Steyn  (1978)]  that  testing  Hq  :  the  population  mean  vector  p^_  of  X, 

X^  ~  ^(p^,  E) ,  remains  constant  during  the  sampling  process  against 
the  alternative  that  p  varies  during  the  process,  is  equivalent  to 


*•  ■>  i  iMhm 
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"k  k 

testing  Hq  :  GV  =  |z|  against  :  GV  =  |Z  |,  where  Z  =  (I  +  2D/n)Z, 

D  =  Zn  Z  \  n  being  the  sample  size.  One  of  the  many  appli¬ 

cations  of  this  result  and  hence  test  for  SGV  can  be  seen  in  multi¬ 
variate  quality  control. 

(b)  Generalized  Canonical  Variable  (GCV)  Analysis.  When  the  original 
vector  can  be  divided  into  k  _>  2  mutually  exclusive  groups,  Anderson 
(1958)  proposed  GCVs  to  be  obtained  by  minimizing  their  GV.  Steel 
(1951)  and  Kettering  (1969)  (in  his  Ph.D.  thesis)  have  attempted  to 
construct  such  GCVs.  However,  no  results  on  statistical  inference 
associated  with  these  GCVs  are  available.  Gnanadesikan  (1977), 
drawing  from  a  well-known  example  in  psychometry,  posed  the  problem 
of  selection  among  GCVs  obtained  by  different  types  of  grouping.  The 
problem  with  its  extension  in  full  generality  boils  down  to  tests  of 
SGVs  and  hence  can  be  tackled  by  the  methods  outlined  in  this  paper. 

(c)  Generalized  Homogeneity  of  Multi-dimensional  Scatter.  Dyer  and 
Keating  (1980)  were  interested  in  the  homogeneity  of  variances  of  the 
sealed  bids  of  'five'  Texas  offshore  oil  and  gas  leases.  A  glance  at 
their  data  reveals  that  the  'five'  leases  are  actually  five  groups 
with  different  number  of  components.  Treating  the  data  as  in  an 
univariate  set-up,  they  proceeded  to  test  for  the  homogeneity  of 
variances.  However,  it  seems  more  appropriate  to  consider  the  groups 
as  vector  variables  and  test  for  homogeneity  of  'variances'  of  t hese 
groups,  which  we  can  term  as  'Generalized  Homogeneity.'  This  will  be 
equivalent  to  testing  homogeneity  of  SGVs,  which  is  precisely 


- 


defined  above. 
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4.  Distributions  and  Associated  Special  Functions.  First  consider 

2pi 

general  Let  s  be  sample  sums  and  product  matrices,  N.  g.^  = 

|s.|,  being  the  sample  size,  i=l,...,k.  Consider  the  setup  of 
Section  2.  Let  be  the  critical  region  of  the  LRTs  or  equivalent 
criteria  for  testing  H  against  both-sided  alternatives,  j=  1,2,3. 
Two-stage  maximum  likelihood  estimation  and  some  judicious  transfor¬ 
mations  are  employed  to  provide  the  following  results. 

l/pi 

:  X^  =  |s_J  >  a^  or  <  a^;  a^,a^  being  constants. 


l/pi  l/p. 

<*>2  :  ^2  =  lSjJ  /lSjl  J  >  b^  or  <  bi*  bo,bl  bein8  constants. 


k  „  „  N.p./2 

W3  :  ^3  =  .n  ^SiO^“NiPi  Si(/SNiPi^  *  *  <  n0  (constant) • 


Through  the  use  of  Calculus  of  Residues,  the  exact  null  and  non-null 
distributions  of  X^  is  given  in  computable  forms  of  G-f unctions , 
those  of  X2  in  computable  forms  of  H-function  and  those  of  u  =  £X^, 

£  being  a  constant,  for  all  and  k  but  p .  = 1  or  2  as. 


fv,a  (u)  =  tr(|)/nr(|  aj)]C^/2(IlaJ)[(2Tr)(n"1)/2/r(^)]uk-1 


(-log  u)  (n  3^2  Ka  (u)  ,  0  <  u  <  1 


where  £  (u)  =  1  +  £  ,  v  (a) (-log  u)r  if  e  2lT^  <  u  <  1  and  v  (a),  a 

a  r— i  r  -  r  —  * 

function  of  Bernoulli  numbers,  v.  C  .  and  6  are  obtainable  when  N  .  k 

a  x* 

and  p^  s  are  specified.  Tables  for  the  distribution  of  u  are  given 

2 

in  Dyer  and  Keating  (1980).  Usual  large-sample  X  approximations  hold 
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for  -2  log  X  ,  j  =  1,2,3.  Additionally,  normal  and  Gamma  approxima- 
C2  3  i/Pi 

tions  for  §£(/^i*  are  rev:i-ewe^*  An  F-approx ima t ion  is  provided 
for  the  distribution  of  X^.  For  X^,  an  approximation  is  given  in  terms 
of  Bartlett's  distribution.  These  approximations  turn  out  to  be  exact 
when  p^ s  are  any  combinations  of  1 s  or  2 s. 

Next  the  equi-correlated  case  is  considered.  The  LRT  for 
HQ  :  pi  =  piQ  against  :  p,^  ^  p^  (i  fixed)  is  derived.  The  MLE  for 
p^  is  shown  to  have  some  undesirable  properties.  A  test  based  on 
truncated  Best  Unbiased  Estimator  for  p^  is  provided.  The  distribu¬ 
tions  of  the  test  statistic  for  both  the  null  and  non-null  cases  are 
represented  in  terms  of  Rummer's  function  and  Whittaker's  function. 

The  LRT  for  is  given.  Additionally,  a  simpler  test  is  provided 
through  a  characterization  of  the  problem  in  terms  of  the  smallest 
eigenvalue  of  the  covariance  matrix.  LRTs  and  modified  tests  using 
Isotonic  regression  techniques  are  also  provided  for  and 

Total  and  partial  unbiasedness  of  most  of  the  tests  discussed 


above  are  also  established. 
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ABSTRACT 


In  1944,  A.  Selberg  [4]  evaluated  an  Important  multivariate  extension  of 
the  beta  integral.  He  showed  that 


(1) 


/  —  /*  fr  r<‘V'l  TT  (trt)|2kdtl  •••  dt 

0  0  U1  1  1  lsi<jsn  1  *  1  n 


_n_  r(x+(J-l)k)  r(y  +  (J-l)  k)  r(l  +  Jk) 
J=»l  r(x  +  y  +  (n+J-2)k)  r(l  +  k) 


where  Re  (x)  >  0  ,  Re  (y)  >0  and  Re  (k)  >  -  i  ,  -  BSizl  . 

n  n-i  n-i 

Upon  normalization,  (1)  gives  the  joint  density  function  of  our  principle 
object  of  study:  the  Selberg  distribution  with  parameters  x,  y,  k  and  n  . 
Some  important  limiting  distributions  are  as  follows.  Let  x  »  y  and  y  —  ® 
Then  (1)  becomes 

n 

2 


(2) 


1 


(2tr) 


a 

f 


--  y  i2 


-00 


"IT  (t,-t 

1  s  i<j  s  n 


2k 

i  ‘j-1  dti  •••<**» 


n  r(l  +  jk) 

3  IT  - 

j  =  l  r{i  +  k) 


Re(k)  >  -  £ 


which  was  studied  by  Mehta  and  Dyson  [3]  .  For 


k  =  j  ,  k  =  1  and  k  =  2  ,  this  corresponds  to  the  distribution  of  the 
eigenvalues  of  orthogonal ,  hermitian  and  symplectlc  matrices,  respectively. 


Letting  y— »  In  (1)  yields 


,oo  oo  n  „  ,  -t. 

/  •••  /  TT  tf'e  1 

0  0  1=1  1 


TT  <t1-t|)l2Ic  dt.  -  -  *  dtn 

lii<jsn  J  1  n 


n  r{x  +  {J-l)k)  Hl  +  Jk) 

*  TT  - 

j  =  i  r(i  +  k) 


where  Re  (x)  >  0  and  Re  (k)  >  -  ~  *  For  k  =  j  ,  this 

corresponds  to  the  distribution  {see  Anderson  [1;  ch.  13])  of  the  eigenvalues 

of  a  random  matrix  from  the  W  is  hart  distribution  with  mean  vector  0  and 
spherical  dispersion  matrix.  , 

For  n=l,  (1),  (2)  and  (3)  reduce  to  the  univariate  beta ,  normal  and 
gamma  or  chi-square  distributions,  respectively.  The  most  fundamental  of 
these  is  the  beta  distribution ,  since  it  has  2  parameters  and  includes  the 
others  as  limiting  cases.  We  are  led  to  study  the  Selberg  distribution  where 
x,  y,  2k  and  n  are  positive  integers.  For  n  =  l,  this  is  the  distribution 

of  the  Xth  order  statistic  from  a  sample  of  x  +  y  -  1  lid  u(0 , 1) 

random  variables.  We  shall  generalize  this  to  Selberg' s  distribution.  Let 
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The  number  of  ways  in  which  this  can  be  done  is  given  by  the  multinomial 
coefficient 

l  M  1 

1.  •  ••  *  l»x-l,  •  •  •  ,  x-1,  y-1,  •  •  • ,  y-1,  2k,  •  •  • ,  2k 

_ Mj _ 

(n) 

[(X-1)!  (y-1) !  ]n  (2k!)  2 

which  we  denote  by  ( ^ )  .  Let  T (x,  y ,  k,  n )  denote  the  number  of  these 
choices  which  also  satisfy  (inequalities  with  sets  hold  for  all  elements) 

(5)  SA  <  mA  <  Lj  ,  lsisn  and  mi  <  <  mj  •  lsi<  )sn  . 

The  condition  (5)  implies  <  ra^  for  1  a  i  <  J  s  n  .  Let  I(x,y,k,n) 
denote  Selberg's  Integral  (1)  .  Then 


(6) 


I(x,  y ,  k,  n)  _ 
n ! 


/  •••  /  TT  tf 1  (i-t.)y-‘  T r  (tr..)d«.--.dt 

<tl<  ■<  tnd  1=1  1  lsicjsn  J  1  1  n 


corresponds  to  the  distribution  of  the  order  statistics  from  Selberg's  distribution. 
We  have 


I  (x,  y ,  k ,  n) 
n ! 


T(x,  y.  k,  n) 


(7) 


so  that  the  problem  of  evaluating  (1)  Is  equivalent  to  the  combinatorial  problem 
of  evaluating  T (x, y ,  k,  n)  .  Let 


u(l)  <  u(2)  <  *  *  *  <  U(M-1)  <  U(M) 

be  the  order  statistics  of  a  sample  of  M  lid  u  (0, 1)  random  variables. 
Choose  one  of  the  objects  counted  by  T  (x,  y ,  k,  n)  at  random,  with  each 
equally  likely  to  be  selected.  Then 

(8)  t*  (tj.  —  .y 

t 

Is  distributed  as  the  order  statistics  of  Selberg's  distribution  (corres ponding  to 
(6)).  Let  a  «  Sn  be  chosen  at  random,  with  each  permutation  equally  likely 
to  occur.  Then 


d(t) =  (u 


has  Selberg's  distribution.  This  fact  and  (7)  both  follow  by  generating  the 
random  variables  In  (6)  from  an  acceptance-rejection  proce;'  rr 
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"The  Zonal  Polynomials  of  Multivariate  Analysis 
as  Special  Functions" 


ABSTRACT 


1.  Introduction 

The  special  functions  which  have  received  the  most 
attention  in  multivariate  analysis  are  the  zonal  polynomials. 
Modulo  the  orthogonal  group,  they  are  roughly  the  eigen¬ 
functions  of  matrix  multiplication.  They  provide  the 
natural  group  theoretic  basis  for  the  space  of  complex¬ 
valued  analytic  functionals  of  an  n  x  n  matrix  which  are 
both  left  and  right  invariant  under  the  orthogonal  group. 

Since  this  space  is  the  home  of  the  multivariate  Gaussian 
density  as  well  as  others,  its  natural  decomposition  is 
necessarily  of  interest. 

The  ease  of  dealing  with  zonal  polynomials  would  be 
enhanced  considerably  if  the  abstract  eigenfunction  property 
were  augmented  with  a  concrete  and  explicit  workable  formula. 
The  desirability  and  difficulty  of  obtaining  such  an  expression 


--  -■ 


i. 


has  been  discussed  on  several  occasions  (James,  1975), 
(MacLaren,  1975) .  It  is  our  purpose  to  give  such  an 
expression  and  note  how  it  arises  from  considering  zonal 
polynomials  as  special  functions  defined  in  the  group 
theoretic  sense  of  zonal  spherical  functions  on  symmetric 
space.  The  discussion  and  results  presented  here  have 
previously  appeared  in  the  author's  Ph.D.  thesis  (Kates, 
1980) .  An  attempt  has  been  made  to  keep  the  exposition 
at  a  fairly  elementary  level. 

2.  Fourier  Analysis 

We  take  the  viewpoint  of  Fourier  analysis  on  groups. 
This  approach  has  not  previously  been  taken  in  the 
statistical  literature.  It  is  nevertheless  a  quite  natural 
and  direct  approach.  We  proceed  partly  by  example. 

Let  G  =  {g1,g2,g3,g^}  be  a  four  element  commutative 
group.  Let  C(G)  =  (f  :  G-  >C}  denote  the  set  of  complex¬ 
valued  functions  on  G.  Each  function  f  is  defined  by  a 
four-tuple  (f (g^) , f (g2) , f (g3) , f (g4) )  showing  that  C(G)  is 
a  four  dimensional  vector  space.  The  operators  T  ,  which 
act  on  C (G)  are  linear  and  are  defined  as 

T  f (g)  =  f (gg' ) . 

They  have  the  effect  of  translating  the  graph  of  f  by  g' . 


The  commutativity  of  G  implies  the  commutativity  of  the  T  . 


B 


The  Tg  are  individually  diagonalizable  commuting  4  by  4 
matrices  and  so  can  be  simultaneously  diagonalized. 

Their  simultaneous  eigenvectors  are  called  zonal  spherical 
functions  and  form  the  natural  basis  for  Fourier  analysis 
on  C (G) . 

If  G  is  non-commutative  the  T^  no  longer  commute  so  we 
cannot  simultaneously  diagonalize  them.  However,  if  we 
are  only  interested  in  a  subspace  V  of  C(G)  then  it  may  be 
true  that  the  T^  operators,  suitably  restricted  to  V,  are 
commutative . 

Now  suppose  that  G  is  not  a  four  element  commutative 
group  but  rather  the  n  by  n  nonsingular  real  matrices 
regarded  as  a  (noncommutative)  group  under  matrix  multi¬ 
plication  and  inversion.  Within  C(G),  the  complex-valued 
functions  on  G,  let  V  be  that  subspace  consisting  of  those 
functions  which  are: 

(i)  expressible  as  a  power  series 
2 

m  the  n  elements  of  their 
argument 

(ii)  orthogonally  bi-invariant, 

namely  f(XXK)  =  f(X),  for  all 
orthogonal  H  and  K. 

Let  P  be  the  projection  of  the  vector  space  C(G)  onto  the 
subspace  V.  Then  the  operators  PT^  acting  on  V  commute  and 


.... 
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their  simultaneous  eigenfunctions  are  called  zonal 
spherical  functions  of  the  pair  of  groups:  real  non¬ 
singulars  /orthogonals .  They  are  also  called  zonal 
polynomials . 


3.  An  Explicit  Formula 

Let  k  =  (k1,k2,  .  ..,  k  )  where  k^  >  k2  £  ...  s  kn  >  0 
are  integers.  Suppose  A^  =  A^(NTXTXN)  is  the  determinant 
of  the  upper  left  i  x  i  submatrix  of  NTXTXN  where  X  is 
a  nonsingular  n  x  n  matrix  and  N  is  an  n  x  n  random  matrix 
whose  entries  are  independent  Gaussian  random  variables, 
each  with  mean  0  and  variance  1.  The  orthogonal  bi¬ 
invariance  property  of  zonal  polynomials  implies  that  they 
depend  on  x  only  through  XTX  so  we  regard  them  as  functions 
of  XTX  in  what  follows.  Our  aim  is  to  show  that 


k.-k_  k_-k,  k  -k  k 

(*)  Zk(XTX)  =  E(A1X  ^  A2  ...A^!  n_X  Ann) 


by  showing  that  (*)  satisfies  the  eigenfunction  property 
and  forms  a  basis  for  V.  Here  E  means  expectation. 

To  evaluate  any  given  zonal  polynomial  multiply  out 
the  integrand  giving  a  polynomial  in  independent  Gaussian 
random  variables,  each  with  mean  0  and  variance  1.  Since 
the  2mt*1  moment  of  such  a  random  variable  is  ( 2m)  1  /  ( 2mml ) 


and  all  odd  mor 


ts  are  0,  substitution  yields  the  polynomial 
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explicitly.  Alternatively,  one  could  Monte  Carlo  the 

function  k  -k  k  -1c  k 

A,1  2  A  2  3  ...  A  n  . 

12  n 

Since  zonal  polynomials  are  eigenfunctions,  their 
normalization  is  unimportant.  Thus  with  change  of 
normalization  (*)  holds  true  when  the  probability 
distribution  of  N  is  replaced  by  the  uniform  probability 
distribution  over  the  orthogonal  group  (also  called  Haar 
probability  measure) .  This  renormalization  is  such  that 
the  zonal  polynomial  equals  1  when  evaluated  at  the 

*  T 

identity.  We  denote  it  C  (X  X) .  In  one  of  these  two 
forms  the  formula  can  be  used  to  prove  many  additional 
facts  -oout  zonal  polynomials  as  well  as  to  reprove  known 
results  in  a  more  direct  fashion. 
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TOPIC:  Infinite  Divisibility 


ABSTRACT 

This  topic  covers  a  set  of  three  papers  on  the  topic  of  infinite 
divisibility.  Infinite  divisibility  has  always  been  tucked 
away  in  the  jargon  of  complex  variables.  If  we  stick  to  integer 
valued  variables,  then  this  tipic  can  be  handled  through  very 
simple  algebra  and  through  the  use  of  computers  whenever  neces- 
ary.  A  necessary  and  sufficient  condition  is  that  a  series  of 
determinants  be  nonnegative.  Thus  the  user  can  test  for  infinite 
divisibility  be  computing  a  bunch  of  these  determinants  and  see 
if  they  are  nonnegative.  In  the  papers,  we  give  proofs  of 
infinite  divisibility  using  these  determinants.  We  also  give 
sufficient  conditions.  One  sufficient  condition  is,  "If  5V  , 
i  0,1,...  are  such  that  P^/P^  ^  are  monotone  nondecreasinc , 
then  the  distribution  described  by  P.  is  infinitelv  divisible." 
A  side  issue  resulting  from  this  extreme  simplicity  and  "computer 
izing  of  the  problem  anti  pulling  it  out  of  the  ranae  of  the 
ibstract  complex  variables"  is  that  we  are  now  trying  to  test 
for  independence  of  observations  in  successive  plots  in  a  fare 
using  the  first  few  determinants  as  our  test  statistic. 
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"The  Communality  Problem  for  Stieltjes  Matrices" 

ABSTRACT 

The  Communality  Problem  is  very  important  in  Factor  Analysis. 

The  problem  is  that  of  reducing  the  diagonal  elements  of  a  given 
correlation  matrix  so  that  the  resulting  matrix  will  be  positive 
semidefinite  and  of  minimum  rank.  The  new  diagonal  elements  are 
called  the  communalities.  We  defind  a  correlation  matrix  as  a 
symmetric  and  positive  semidefinite  matrix  with  unities  along 
the  diagonal  and  fractional  numbers  between  -1  and  +1  in  the  off- 
diagonal  positions. 

Many  researchers  have  studied  the  Communality  Problem  (for  a 
detailed  set  of  references,  see  Harman  [ 1 ] ) -  However,  no  effective 
solution  procedures  have  yet  been  devised.  In  this  paper,  we  propose 
a  variant  problem  and  give  an  algorithm  for  its  solution.  We  prove  that 
a  solution  to  our  problem  also  solves  the  Communality  problem  if  the 
given  matrix  is  Stieltjes. 

We  can  state  the  Communality  Problem  as  follows: 

Problem  1^ 

Given  a  correlation  matrix  R,  find  a  non-negative  diagonal  matrix 
D  such  that 

(i)  R  -  D  is  positive  semidefinite,  and 
(ii)  rank  (R  -  D)  =  min.  □ 


_ 
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Let  us  examine  the  variant  problem: 

Problem  2 

Given  a  correlation  matrix  R,  find  a  non-negative  diagonal  matrix 
D  such  that 

(i)  R  -  D  is  positive  semidef inite, 

(ii)  rank  (R  -  D)  =  min.,  and 

(iii)  rank  (R  -  D)  +  rank  (D)  =  rank  (R) .  □ 

Condition  (iii)  is  the  additional  constraint  that  we  have  imposed. 

For  convenience,  we  use  the  following  notations: 

Notations 

We  let  CP (A)  denote  Problem  1  with  A  as  the  given  correlation 
matrix,  and  let  P(B)  denote  Problem  2  with  B  as  the  given  correlation 
matrix.  □ 

For  a  given  correlation  matrix  R,  let  and  be  the  nonnegative 
diagonal  matrices  that  solve  CP(R)  and  P(R),  respectively .  The  extra 
condition  for  P(R)  implies  that 

rank  (R  -  D^)  <  rank  (R  -  D2) . 

We  can  easily  construct  an  example  for  which  the  inequality  above  is 
strict . 

Let  us  describe  an  algorithm  for  solving  P(R),  where  R  is  a  given 
correlation  matrix.  Let  an  eigenvalue  decomposition  of  R  be 

R  =  qzqt. 


A 


where  Q  is  orthogonal  and 


(:■:) 


with  a  nonsingular  diagonal  matrix  of  order  r. 


Q  =  (Qr  Q2), 


where  is  n  x  r. 
Define  the  index  set 


Z  =  l  {i  |  row  i  of  Q.  consists  of  all  zeros),  if  r  <  n. 


(1,2, .. . ,n) 


if  r  =  n. 


If  Z  is  empty,  we  cannot  reduce  the  rank  of  R  without  giving  the 
resulting  matrix  a  negative  eigenvalue.  Hence  we  assume  that  Z  is 


nonempty.  Let 


Z  -  U(l).  i (2) . i(£)}. 


For  j  “  1,2,  ....  £,  we  can  find  a  vector  ^  such  that 

^KJ)  '  ~i(j) ' 

The  general  solution  is 

Si(j)  "  Q1E1  Ql~i ( j )  + 


where  w  c  N(R) . 

Let  Y  be  an  d  x  d  matrix  with  its  (j,k)  element  equal  to  the  i(J)-th 
element  of  Note  that  the  elements  of  Y  are  not  affected  by  the 

choice  of  w  in  the  previous  equation,  as  w  .  .  =  0  for  j  «  !,...,£. 


Let  I  be  the  set  of  indices  such  that  for  i,j  e  I 


0 


for  i  *  j. 


Construct  the  diagonal  matrix 
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jii°i(j)Si(j)~i(j)* 


where 


Note  that  is  positive  because 

y.  -  e  Q  Z  “1Qte 

U  ~i  0)1  1 


and 

%(j)  1  N(R)* 


it  Is  easy  to  check  that 


(R  -  D)x1(J)  *  S  for  j  e  1. 

We  can  also  show  that  the  matrix  R  -  D  is  positive  semidefinite  and  that 

rank  (R  -  D)  *  rank  (R)  -  rank  (D). 

As  our  goal  is  to  minimize  rank  (R  -  D) ,  we  want  to  determine  the 
index  set  I  of  maximum  size.  It  is  known,  however,  that  determining 
such  an  I  requires  work  which  is  exponential  in  1,  the  size  of  Z. 
Fortunately,  1  is  usually  very  small  for  the  matrices  of  the  Com- 
munallty  Problem. 

We  next  consider  the  case  when  our  correlation  matrix  is  a 
Stieltjes  matrix,  l.e.  a  positive  definite  matrix  with  nonpositive 
off-diagonal  elements.  Let  D  be  a  nonnegative  diagonal  matrix  which 
solves  P(R).  We  are  going  to  show  that  D  is  also  a  solution  to  CP(R). 
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"The  Laguerre  Transform" 


ABSTRACT 

A  novel  transform  is  presented  which  maps  continuum 
functions  (such  as  probability  distributions)  into  disc¬ 
rete  sequences  and  permits  rapid  numerical  calculation 
of  convolutions,  multiple  convolutions,  and  Neumann 
expansions  for  Volterra  integral  equations.  The  trans¬ 
form  is  based  on  the  Laguerre  polynomials,  associated 
Laguerre  functions,  and  their  convolution  properties. 

Part  1  of  this  paper  deals  with  functions  having 
support  only  on  [0,°°)  .  The  resulting  unilateral 
Laguerre  transform  finds  applications  in  convolution  of 
such  functions,  inversion  of  Laplace  transform,  and  in 
solution  to  renewal  and  related  Volterra  integral 
equations . 

Part  2  of  this  paper  deals  with  functions  having 
support  on  (-<*,«>)  via  a  bilateral  Laguerre  transform 
which  is  an  extension  of  the  unilateral  transform. 
Applications  of  this  technique  include  convolution  of 
such  functions  and  analysis  of  the  Lindley  process. 

Part  1  has  been  published  in  Applied  Mathematics 
and  Computation  and  part  2  has  been  submitted  for  pub¬ 
lication  in  that  journal. 


SUMMARY  “  PART  1 


I 


One  often  encounters  in  applied  studies  integral  equations  [16] 
either  of  form 


a (x  -  x' ) f  (x1 )dx' 


»  b(x) 


or  of  the  form 


(1) 


f  (x) 


a(x  -  x')f{x')dx' 


b  (x) 


(2) 


where  a{x)  and  b(x)  are  specified  functions  and  f(x)  is  to  be 
found.  Equations  (1)  and  (2)  are  said  to  be  Volterra  integral  eq¬ 
uations  cf  convolution  type  cf  the  first  and  second  kind  respective¬ 
ly.  The  Neumann  series  solution  of  (2)  has  the  form.  [19! 


or 


i 


(3) 


where  the  asterisk  denotes  convolution  art 


is  the  k 


Id 


convolution  of  a(x' 

cc 

The  entity  X 
C 

systems  of  integral 
tions  research  [  c  1 , 


with  itself. 

ty ) 

a  *  (x)  and  matrix  variants  associated  with 


equations  of 
encmeermr 


convolution  type  arise  in  oper?- 
r  •.  "i ,  end  biological  studies  T *  n . 
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Sometimes,  differential-integral  equations  give  rise  to  expressions 
such  as 


(At)* 

(k  +  1)  1 


a(k+ll(T) 


which  describes  the  busy-period  density  for  certain  M | G | 1  queueing 
systems  [18]. 

In  easy  cases  the  integral  equations  may  be  solved  analy¬ 
tically  via  Laplace  transformation,  and  full  answers  may  be  obtained 
when  the  Laplace  transforms  are  invertible.  More  often  than  not, 
such  transforms  cannot  be  inverted  and  expressions  such  as  (4)  are 
of  limited  value  when  they  cannot  be  evaluated  explictly.  The 
Laguerre  transformation  techniques  developed  in  this  paper  may 
then  be  of  value. 

The  deconvolution  problem  of  finding  f(x)  from  (1)  when 
a(x)  ana  b(x)  are  known  numerically,  say.  is  particularly 
troublesome,  and  start-up  difficulties  described  belov  may  make 
conventional  numerical  procedures  useless. 

The  Laguerre  transform  techniques  described  map  continuum, 
functions  into  sequences,  and  map  the  continuum  convolution  opera¬ 
tion  into  lattice  convolution  of  these  sequences.  Such  discrete 
convolutions  are  well  matched  to  modern  computer  competence,  and 
mapping  back  to  the  continuum  is  direct. 


the  inversion 
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Laguerre  transformation  has  been  developed  as  a  tool  for  the 
solution  of  differential  equations  [12]*  The  applications  of 
interest  here  are  quite  different  and  new  tools  have  been  needed  to 
convert  the  underlying  simple  idea  into  a  flexible  working  proce¬ 
dure  adapted  to  computer  requirements. 

The  first  section  introduces  the  Laguerre  transform 

T  -f  CD 

T:  f(t)  -*•  (f^  )Q  in  a  form  convenient  for  our  needs.  One  has 

CO 

f  CT J  *  2Z  (5) 

n=0  n  n 

for  any  square-integrafcle  function  f (~ )  on  (0,“',  where  = 

-t/2 

Ln(T)e  are  the  classical  orthonormal  Laguerre  functions  and 

Ln(t)  are  the  Laguerre  polynomials.  The  notation  of  Abramovitz 
and  Stegun  [  1  ]  is  employed  throughout.  Or th enema  1  d  ty  provides 
the  inverse  transformation 


f  f(Tj£  (T)dT  • 


:e) 
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where  $(s)  is  the  Laplace  transform  of  f(T>.  This  relationship 
permits  evaluation  of  f^+  for  many  important  f(t). 

Section  2  provides  simple  examples  of  the  transform,  and 

Section  3  discusses  the  structure  of  T^(u)  in  the  complex  u- 

plane.  Such  insight  into  structure  in  the  complex  plane  is  crucial 
to  many  of  our  algorithms  and  theorems. 

Algorithms  for  the  calculation  of  the  Laguerre  coefficients 
are  presented  in  Section  7.  Section  8  is  devoted  to  a  discussion 
of  the  deconvolution  problem. 

A  variety  of  numerical  examples  of  the  method  are  treated  ir. 
Section  S,  and  the  implementation  of  the  procedure  is  discussed. 

Section  10  describes  interpolation  methods  and  problems  when 
the  known  functions  are  known  only  numerically. 

A  final  section  deals  with  possible  generalisations  c:  the 
method  tc  special  families  of  functions. 
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SUMMARY  -  PART  2 

In  a  previous  paper,  hereafter  designated  by  [A],  a  description  was 

given  of  a  Laguerre  transformation  which  maps  a  function  f(t)  in  L2(0,<=) 

into  a  sequence  (fn)g  on  the  nonnegative  integers.  Moreover,  for  two 

such  functions  f(t),  g(t),  the  convolution  f(t)g(T)  is  mapped  into  the 

#  ^  #  #  ft 

lattice  convolution  (f*g)n  =  I  *n-m^m  =  Cne  °btains  thereby 

®  (kl 

an  algorithmic  basis  for  the  computation  of  multiple  convolutions  fv'J(-r) 

and  related  infinite  series  of  importance  to  statistics  and  applied  proba¬ 
bility. 

Such  Laguerre  transforms  have  one-sided  functions  as  their  natural 

domain  because  the  Laguerre  polynomials  Ln(t)  and  Laguerre  functions 
-  /  *■> 

l  (t)  =  L  (T)e’  *"  are  associated  with  the  one-sided  weight  function 

e”T  on  (O,™).  Nevertheless,  the  methods  have  a  simple  extension  to  two- 
sided  functions  on  the  full  continuum  (-■*,  «)  via  the  same  Laguerre 
functions  as  we  will  see. 

A  variety  of  applications  exist  to  statistics,  operations  research, 
and  engineering.  In  statistics,  for  example,  one  has  need  for  multiple 
convolutions  cf  two-sided  distributions  unavailable  analytically ,  that 
of  the  logistic  distribution,  for  example.  Even  relatively  innocuous 
distributions  such  as  the  Laplace  distribution  convolve  with  difficulty. 

In  operations  research  studies  dealing  with  queues,  inventories  ar.d 
storage  systems,  one  encounters  as  a  structural  entity  [5]  the  extended 

ejo 

renewal  densitv  h(x)  =  )  a^(x),  where  afx)  is  a  probabilitv  density 

1 

function  with  two-sided  support.  Fcr  many  densities  of  interest,  evalua¬ 
tion  of  h(x)  has  beer,  resistant. 


B-79 


In  the  earlier  paper  [A],  the  crucial  role  of  the  complex  plane  in 
the  formulation  of  the  algorithms  was  evident,  even  though  the  algorithms 
were  entirely  in  the  real  domain.  For  the  bilateral  transform,  the  com¬ 
plex  plane  is  again  very  much  present,  with  Laurent  expansions,  bilateral 
Laplace  transformation  and  conformal  mapping  entering  as  crucial  tools. 

The  first  section  extends  the  earlier  formalism  to  the  full  continuum. 
That  this  extension  is  natural,  and  not  just  an  artificial  piecing  together 
of  the  formalism  for  each  half-line,  will  be  clear  from  (1.9),  (1.12)  and 
(1.13).  The  harmony  of  the  basis  will  also  emerge  vividly  in  Section  3, 
which  deals  with  the  extent  of  the  transform  coefficients,  and  associated 
uncertainty  relations.  The  topic  of  extent  is  crucial  to  the  utility  of 
the  Laguerre  transform  method  as  a  numerical  tool.  Numerical  examples  are 
presented  in  Section  5.  A  table  of  contents  provides  the  reader  with  an 
overview  of  the  paper. ' 


Two  references  (V.  I.  Krylov  and  N.  S.  Skoblya  [8],  and  K.  T.  Keeks  [12]) 
have  come  to  the  authors'  attention  subsequent  to  publication  of  [A],  Beth 
deal  with  the  use  of  Laguerre  functions  for  the  numerical  inversion  of  c-r 
sided  Laplace  transforms. 


Paul  Speckman 
University  of  Oregon 


"Spline  Smoothing  and  Optimal  Rates  of  Convergence 
in  Nonparametric  Regression  Models" 

ABSTRACT 

1.  The  model.  Consider  the  nonparametric  regression 

model 

yi  =  f(xi.)  +  ei  »  i  *  1, .  •  .  ,n, 

where  observations  are  taken  at  distinct  points  assumed  for 
simplicity  to  be  in  [0,1].  The  usual  assumptions  on  the 

2 

random  errors  are  in  force,  i.e.,  Ee.  =  0,  Ec.e.  =  j  °  > 

but  the  response  function  f  is  assumed  only  to  be  sufficientl 

smooth  so  that  llf^^lt^  ».  J  f^^(x)^dx  exists  and  is  finite. 

0 

This  model  is  motivated  by  certain  robustness  consider¬ 
ations.  For  small  a  >  0,  the  class  if:  f  has  k-1  abs.  cont. 
derivatives,  life'll  s  a]  can  be  viewed  as  a  collection  of 
response  functions  at  least  locally  well  approximated  by 
polynomials  of  degree  k-1  (or  order  k) .  If  a  regression 
method  is  uniformly  good  in  this  class,  it  is  robust  to 
arbitrary  small  departures  from  the  standard  kth  order 
polynomial  model  (see  [7]).  This  concept  is  also  related  to 
the  models  of  Sacks  and  Ylvisaker  [6]. 

2.  A  basis  for  splines  and  the  proposed  estimator.  A 
variant  of  spline  smoothing  is  proposed  which  is  most  easily 


I 

A 


expressed  using  the  basis  of  Demmler  and  Reinsch  (2].  Given 
a  distinct  knot  set  (x^ , . . . ,xn) ,  let  denote  the  n  dimen¬ 

sional  space  of  natural  polynomial  splines  of  degree  2k-l 
wit.i  simple  knots  at  the  prescribed  points.  There  is  a  basis 
. s  j  1  for  along  with  eigenvalues  t  X ^  j  ^  =  ^ 

determined  (essentially  uniquely)  by  the  conditions 

n 


cp .  (x.  )cp .  (x,  )  =  f .  . 
t  f  VJ  ?  i  J 

I  v-  (x)cp .  (x)d::  =  '..a. 

”  a  V  /VJ  V  1J  J 


)  -J 


.ii. 


1a  i  0.  -  £  y.a).(x.).  (Note  that  h.  is  the  least  squares 

1  i=L  ^  3  L  3 

n  2  k 

e s i  i motor  of  j%.  if  f  =  f-CP*  €  >  .  )  With  this  basis. 

J  j-i  11 

the  f ..mi J y  of  estimators  to  be  studied,  indexed  by  a  parameter 
X  0,  has  the  form 

n 


:x(x)  -  E^(l-VxXi)+tsicpi  (x)  . 


One  can  see  that  f^  is  a  natural  polynomial  spline,  but  it 
is  not  the  smoothing  spline  of  Reinsch  [5]. 


3.  Optimal  rates  of  convergence.  The  purpose  of 

A 

introducing  f^  is  that  it  is  minimax  in  the  following  sense. 

A 

Let  Ll  be  the  class  of  all  estimators  f  which  are  linear 

A 

in  the  observations.  (Clearly  f.  €  C.)  For  fixed  f  and 

n  * 

f  C  define  T(f,f)  =  —  ■£  (f  (x.  )-f  (x.  ))^.  Then  the  next 

n  i=l 

result  follows  from  Kuks  and  Olman  [4] . 
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A  * 

THEOREM  1.  min  max  ET(f,f)  =  min  max  ET(f,f^) 

f  €  C.  |if(k)l!  *  a  X  *  a 

=  min  £a2\  +  o2  £  (1-,/U.)2}  . 

X  i=l  1  + 

From  the  last  expression,  the  minimax  rate  of  convergence 
of  ET(f,f^)  is  seen  to  depend  critically  on  the  eigenvalues 
.  )i  j  .  In  this  paper  an  approximation  is  given  when  the 

x . ' s  are  equally  spaced,  although  similar  results  undoubtedly 
ho i u  i or  any  setting  where  the  points  are  suitably  regular. 

■  i.>  kt  y  tool  is  an  estimate  first  given  by  Utreras  [  -  j  and 
due  iu  the  form  here  to  the  author. 

LEMMA.  For  k  >-  2,  there  exist  constants  d,  such  chat 
_2k  K 

X j  -  d^  (1  +  o(l))  uniformly  in  j  for  j  =  0(n2//  ^2k  +  ^ ) . 

With  this  estimate  the  best  possible  uniform  rate  of  convergence 

follows. 

THEOREM  2.  For  equally  spaced  points  and  k  &  2, 

min  max  ET(f,f)  -  n"2k/ (2k  +  X)c,  (a,  o)  (1  + o(l)> 

f€C  !lf<k)Ma 

where  ck(a,a)  is  a  constant  depending  only  on  k,  a,  and  o. 

This  rate  compares  with  recent  work  of  Stone  [8].  There 
the  same  rate  is  obtained  under  weaker  conditions  but  without 
the  exact  constant. 

In  the  usual  situation  there  is  a  single  fixed  but 

*  A 
unknown  f.  Let  X  be  the  minimizer  of  ET(f,f^). 


LEMMA.  ET(f,f  *)  -  „-2k/(2k+l)  c  (||f  (k>||  >0)  (l  +  Q(l))  and 
X  K 

X*  =  o^1^2^1)). 

4.  An  adaptive  estimator.  Since  in  practice  neither 
! i  f  ^ k ^ ;  nor  o  is  generally  known,  an  experimentor  is  forced 
to  choose,  subjectively  or  otherwise,  some  value  of  the 
smoothing  parameter  X  to  use.  This  problem  is  characteristic 
of  virtually  every  nonparametr ic  regression  method.  The  pro- 
i  lu  re  is  to  estimate  X  by  Wahba' s  method  of  genera  1  i /.t. 

crass  validation.  The  main  new  result  is  weak  consistency  and 
asymptotic  optimality  of  the  resulting  adaptive  method.  This 
procedure  is  closely  related  to  the  one  in  Craven  and  Wahba  []] 
where  ordinary  smoothing  splines  are  used.  However  the  con¬ 
sistency  result  here  is  much  stronger. 

a 

Since  is  a  linear  estimate,  there  is  an  n  x  n 

A  A  A 

matrix  A(X)  such  that  =  (f^(x^)  ,  .  .  . ,  f^(xn))  '  =  A(X)y. 

Then  the  cross  validation  function  is  defined  to  be 

VX)  *=  H,,y'fXl|2/I  nTr  (I"A(X))]2  • 

Also,  redefine 

Vl>  =T<£"V- 

Then  a  straight  forward  application  of  the  GCV  theorem  of 
Golub,  Heath,  and  Wahba  [3]  yields 


uniformly  in  0  £  X  £  X^  ,  where  X°  is  some  sequence 
^  o 

satisfying  Xn/Xn  0*  This  implies  that  the  minimizers  of 
EV^(X)  and  ETn(X)  asymptotically  coincide.  If  in  addition 
normality  is  assumed,  the  following  much  stronger  result  is 
true . 


LEMMA .  Suppose  the  distribution  of  y .  is  normal  for 
i  -  l,...,n.  Then  there  is  an  estimator  nn  of  c~  such 
i  ha  t 


Vn(X)-o/-ETn(»  _ 

ET  (X)  °PU) 

n 


uniformly  in  0  ^  X  ^  X°  .  From  this  the  main  result  follows, 
- n 


THEOREM  3.  Let  X_  be  the  minimizer  of  ET_(X).  Under  the 
-  n  -  n  - 

assumptions  above,  there  is  a  sequence  of  (possibly  local) 

*  -  *  P 

minimizers  { X_3  of  V_(X)  such  that  X_/X  -»  1,  Moreover, 

-  n  d —  n  -  n  n  - 


Tn<Xn>/EW 


1. 
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"Incomplete  Gamma  Functions,  Numerical  and  Asymptotical 
Aspects  for  Evaluation  and  Inversion" 


ABSTRACT 


following  topics  will  be  discussed: 

-  The  asymptotic  expansion  of  integrals  of  the  type 


: : ) 


F.  (X) 


(a/t } 


f  (t) 


dt. 


a  -*■  °° 


•  'nich  is  uniformly  valid  with  respect  to  x  in  a  domain  that  contains 

x  =  0. 

-  The  inversion  (for  large  a)  of  F  (x) ;  i.e.,  the  computation  of  x  from 
the  equation  F  (x)  =  q,  where  q  and  a  are  given  (a  large). 

cl 

-  The-  numerical  evaluation  of  the  incomplete  gamma  functions  for  large 
values  of  their  parameters.  The  algorithm  is  based  on  the  asymptotic 
representation  of  F  (x) . 

-  Tne  inversion  of  the  incomplete  gamma  functions  for  large  values  of 
tneir  parameters. 


We  suppose  that  in  (1)  f  is  analytic  in  a  domain  of  the  complex 
t-plane  that  contains  the  real  axis,  and  that  f ( t)  >0  for  real  (finite) 
t.  Furthermore  we  suppose  that  F  (-<*>)  =  1.  Several  distribution  function 
of  mathematical  statistics  can  be  written  as  (1),  for  instance  the 
incomplete  gamma  functions.  The  function  f  may  depend  on  a;  we  are 
allowing  a  representation 
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®  f  (t) 

f(t)  ~  £  n  -  ■  a  -*•  ®, 

n=0  a 

in  which  f  and  f  have  the  same  domain  of  analyticity. 

The  asymptotic  expansion  of  F  (x)  is  obtained  by  an  integration  by 

cl 

parts  procedure.  It  has  the  form 


(2) 


F  (x> 

a 


erfc (xa^) 


00 


I 

n=0 


2 

-ax 

e 


00 


I 

n=0 


Bn(x) 

n+1  ‘ 

a 


The  A  and  B  are  defined  recursively  in  terms  of  f.  ^ey  have  the  same 
n  n 

domain  of  analyticity  as  f. 

The  inversion  of  Fa(x)  =  q,  0  <  q  <  1,  is  for  large  a  based  on  (2) . 
First  we  solve  by  using  known  algorithms  the  equation 

l 

2  erfc(xa-)  =  q, 


giving  for  x  an  approximation  x^.  Then  an  asymptotic  expansion 

x  ~  xo +  K  +  "1*2 + 

a 

is  derived,  in  which  the  are  expressed  in  terms  of  f,  for  instance 


xi  = 


In » 1  +  f(xQ)  -  f (0) ) 


i 


The  expansion  is  especially  useful  in  the  neighborhood  of  x  =  0,  i.e., 

for  valuer  of  q  near  Numerical  experiments  for  the  case  of  the 

incomplete  gamma  functions  show  uniformity  of  (3)  with  respect  to 

q  e  [0,1].  For  these  functions  more  information  on  x  will  be  given. 

n 

In  general,  the  expansion  in  (2)  is  too  complicated  for  numerical 
computations .  This  will  be  illustrated  for  the  incomplete  gamma  func- 


for  a  and  x  large  near  the  critical  line  x  =  a.  The  algorithm  runs  as 
follows.  Introduce  for  x  i  0  and  a  >  0 


B-88 


(5)  X  =  x/a,  n  =  {2[X-l-ln  X]}^ 

with  sign  n  =  sign(X-l)  (X  >  0).  Then  (see  [1]) 

P(a,x)  =  \  erfc[-  (a/2)  ]  -  R  (n) 

a. 

(6) 

Q  (a ,x)  =  \  erfc[  (a/2)^n]  +  R=(n) 

a 

and  R  (n)  has  an  expansion  as  the  second  series  in  (2).  The  function 

a 

S  (n)  in 

a 

:  _i  2 

R  in)  =  (a/2-)2  e  2an  S  I n) 
a  a 

is  for  a  >  0  and  -•  1R  slowly  varying  and  it  is  analytic  in  a  neighborhood 
of  ti  =  0 : 


-7) 


JC 

s  (n)  =  ^  s  (a) nk 

a  ,  _  K 

K  —0 


Ini  <  2  r 2  . 


It  satisfies  the  differential  equation 


It  Vnl 


+  a  s  Cn) 

a 


r*(a) 


n 

X-l 


★  .  p  cl  —<3. 

where  T  (a)  =  (a/2Ti)  c  e  a  T(a)  (a  >  0)  ,  the  relation  between  n  and  \ 

being  given  in  (5).  The  s^  in  (7)  are  easily  obtained  from  a  recursion 

relation,  which  for  numerical  applications  is  used  in  backward  direction. 

Instead  of  (7)  we  can  expand  S  (n)  in  Chebyshev  polynomials 

a 

T^ (n/ng) »  for  some  Hg  >  0,  yielding  an  even  beter  expansion.  In  both 

cases  we  obtain  expansions,  which  converge  faster  as  a  increases,  and 

from  which  S  (n)  can  be  computed  for,  say,  -1  <  n  <  1,  or  equivalently 
a 

for  0.3017..  <  X  <  2.357  ...  .  This  gives  an  algorithm  for  the  functions 
in  (4)  for 


0.3017a  <  x  <  2.357a. 


Initially  we  supposed  large  values  of  a.  The  algorithm  works  quite  well 
for  a  £  5. 

[1.1  N.M.  TEMME,  The  asymptotic  expansion  of  the  incomplete  gamma  functions , 
SIAM  J.  Math.  Anal.  10,  1979,  757-766. 


«*• '  i.  .. 
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"Applications  of  Sphor  '  if:  ; Wave  Functions  to 
Time  Series  Ana! vs  is" 


ABSTRACT 


(iiven  a  Unite  sample  of  a  process  some  of  the  major  problems  of  time 
series  analysis  are  those  of  testing  for  the  presence  of  a  line  component,  choos¬ 
ing  an  algorithm  to  estimate  the  spectrum  so  that  the  estimate  is  not  dominated 
by  bias,  ensuring  that  the  estimate  is  consistent  and  statistically  meaningful, 
and  maintaining  these  properties  in  the  presence  of  minor  variations  of  assump¬ 
tions.  Despite  a  long  history  these  problems  are  still  lacking  satisfactory  solu¬ 
tions  in  till  but  the  simplest  cases. 


We  assume  a  finite  sample  { jt0 ,  x,  ,  •  •  •  ..viV - 1 }  of  a  wide  sense  station¬ 
ary  time  series  having  the  centered  Cramer  representation 

X„  =?  J  C  dZ  (  v) 


The  extended  Cramer  representation  permits  a  distinction  between  harmonic 
analysis  and  spectrum  estimation:  harmonic  analysis  is  concerned  with  the  firs! 
moments  of  <//  (i<),  while  spectrum  analysis  is  the  problem  of  estimating  the 
second  moments  of  dZ  (  r ). 

These  moments  arc  estimated  as  functions  of  the  discrete  Fourier 
transform  of  the  observations  which,  for  notational  simplicity,  it  is  convenient 
to  define  in  centered  form: 


tv  a  -i  2  »/(/i- 

X{f)  =-■  2/  C  *a 

n^O 


Using  the  spectral  representation  for  the  data  in  this  formula  we  have 


V; 


*(/>“  / 

-  b 


sin  .V  7t  ( /  —  i’  ) 
sin  7r  (/  —  i’ ) 


dZ(r) 


(1) 


which  is  the  convolution  of  the  Cramer  process,  dZ  (  r ),  with  a  Dirichlet  ker¬ 
nel. 


In  these  paper  we  give  a  new  solution  to  these  problems  obtained  by 
applying  a  ''localized”  Karhunen-Locve,  or  principle  components,  expansion  in 
the  frequency  domain  to  estimate  the  moments  of  the  Cramer  process,  dZ  (  r). 
From  this  viewpoint  equation  (1)  is  best  regarded  as  a  linear  Fredholm  integral 
equation  of  the  first  kind  for  dZ  ( v )  and,  since  detailed  information  about  the 
eigenfunctions  and  eigenvalues  of  the  Dirichlet  kernel  have  recently  been  pub¬ 
lished  by  Slcpian  1 1 978],  it  is  feasible  to  attempt  it’s  solution.  These  eigenfunc¬ 
tions,  denoted  by  Uk  (N.  W\f  ),  k  =0.  1,  •  •  •  A'  — 1  arc  known  as  discrete  pro¬ 
late  spheroidal  wave  functions  and  are  solutions  of  the  equation: 
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IV 

f  sll?--//-=^;1  =  W)-Uk{N.  ^l/) 

^  smir(/-/  ) 

Because  the  Dirichlet  kernel  is  degenerate  it  is  impossible  to  obtain  exact  or 
unique  solutions;  what  we  attempt  is  an  approximate  solution  which  is  both 
numerically  and  statistically  plausible. 

Within  this  framework  the  harmonic  analysis  line  test  procedure  becomes 
essentially  an  analysis  of  variance  applied  to  the  coefficients  of  the  cigenexpan- 
sion  and  so  results  in  an  approximate  likelihood  test.  Similarly,  from  the  spec¬ 
trum  estimation  viewpoint,  the  technique  used  to  approximately  solve  the  fun¬ 
damental  integral  equation  results  in  an  estimate  which  is  data  adaptive  and 
computationally  equivalent  to  using  the  weighted  average  of  a  series  of  direct 
spectrum  estimates  made  with  orthogonal  data  windows  (  discrete  prolate 
spheroidal  sequences  )  applied  in  the  time  domain.  Since  the  expansion  is 
applied  in  the  frequency  domain  it  is  insensitive  to  minor  departures  from  nor¬ 
mality  and  the  time  domain  aspects  of  the  procedure  permit  it  to  be  easily 
robustificd  against  gross  outliers. 

While  this  procedure  is  philosophically  very  different  from  the  various 
autoregressive  and  maximum-entropy  methods  currently  fashionable  the 
analysis  of  variance  procedure  provides  considerable  insight  into  the  “super- 
resolution”  question.  Further,  in  addition  to  providing  estimates  of  the  spec¬ 
trum  which  are  based  on  well-established  principles  instead  of  heuristics  this 
methodology  also  permits  a  resolution  of  the  differences  between  windowed  and 
unwindowed  philosophies. 


Slcpian,  D.  [1978]  Prolate  Spheroidal  Wave  Functions.  Fourier  Analysis,  and  Uncertainty-  V:  The 
Discrete  Case,  Bell  System  Tech.  J.  57  pp  1371-1429. 
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1.  Introduction 


ABSTRACT 


In  this  paper,  we  propose  the  use  of  MACSYMA  to  aid  in 
obtaining  absolute  error  bounds  for  a  class  of  asymptotic  ex¬ 
pansions.  The  uniform  asymptotic  expansions  of  interest  here 
are  those  referred  to  in  statistical  literature  as  Edgeworth 
series.  The  Edgeworth  series  can  produce  quickly  converging 
and  accurate  approximations  to  statistical  distribution  func¬ 
tions  whose  computations  arc  ordinarily  intractable.  This  is 
especially  true  for  noncentral  distributions  with  large  para¬ 
meters  which  include  the  noncentral  beta  and  gamma  distribu¬ 
tions.  Unfortunately,  the  only  way  to  determine  the  accuracy 
of  these  asymptotic  expansions  is  by  comparison  to  known 
values . 

The  plan  of  attack  in  determining  error  bounds  is  to  use 
MACSYMA  to  find  a  continued  fraction  ' corresponding '  to  the 
asymptotic  series.  Bounding  an  asymptotic  series  or  summing 
a  divergent  series  via  continued  fractions  is  certainly  not 
new  (Wall  [10],  Henrici  [6],  Shenton  and  Bowman  [9]).  How¬ 
ever,  this  will  be  the  first  time  the  technique  has  been 
applied  to  the  rather  complex  Edgeworth  series.  This  fact  is 
not  surprising  because  the  proposed  technique  would  be  more 
than  formidable  without  MACSYMA. 

2 .  The  Edgeworth  Series 

Various  derivations  of  the  Edgeworth  series  can  be  found 
in  the  statistical  literature.  See  Berry  [1],  Draper  and 
Tierney  [3],  Esseen  [4],  and  Hsu  [8].  We  prefer  an  exposition 
similar  to  that  of  Hill  and  Davis  [7],  because  the  series  is 
presented  in  an  explicit,  easily  programmable  form,  which 
eliminates  the  bother  and  possible  mistakes  in  the  production 
and  use  of  numerous  tabled  constants.  See  Draper  and  Tierney 
[3]. 


Given  a  random  variable,  X,  with  pdf,  f (x j 0) ;  where  0  is  a 

possibly  vector  parameter.  The  characteristic  function  (cf) 
(or  Foiirier-Stielt jes  transform)  of  f(x|0)  is  denoted: 

<f>j,(t|0)  =  E[exp(itx)]  =  /  _ooexp  (itx)  f  (x  |  0)  dx .  (2.1) 

If  a  power  scries  expansion  in  (it)  exists  for 
In  <f>x(t|G)  ,  i.e.  , 
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ln*x(t|0)  »  c(j)(it)j/j!f  (2.2) 

which  is  usually  valid  only  if  |t|  <  1,  then  the  coeffi¬ 
cients,  K^(j),  of  (it)-Vj!  are  the  "cumulants"  of  the  distri¬ 
bution  of  the  random  variable,  X. 

This  starting  point  for  developing  the  Edgeworth  expan¬ 
sion  is  closely  related  to  the  central  limit  theorem,  in  that 
we  can  view  the  random  variable,  X,  as  being  the  sum  of 
2 

n  =  o  =  Kx(2)  independently  and  identically  distributed 
(iid)  random  variables.  The  cf  of  the  sum  of  n  iid  random 

variables  has  the  form  [<{>(t)]n.  Clearly,  if  the  cumulant 

9  2 

expansion  (2.2)  exists,  4>x(t|0)  =  [exp{  (lnij>x  (t  |  0)  ) /a  }  ]° 
which  implies  that  the  random  variable,  X,  can  be  viewed  as 
the  sum  of  n  =  a  iid  random  variables  with  cumulants  Kx(j  ) /a  . 

If  we  make  the  transformation,  Y  =  (X-y)/o,  and  note 
that  the  cumulants  of  the  standardized  random  variable,  Y, 
are : 

Ky(l)  =  Kx(l)  -  y  =  0 , 

ky(2)  =  Kx(2)/o2  -  1,  and  (2.3) 

Ky(j)  =  Kx(j)/o^,  for  j  :>  3; 

we  can  write  the  cf  of  Y  in  terms  of  the  cumulants  of  X: 
*Y(t|0)  =  exp  j-t2/2  +  f3Kx(j  )  (it/a)  j  j  I 

=  *z(t)  •  expj  IiX(j)(it)j+2/(aj(j+2)!)|  ;  (2.4) 

2 

where  X(j)  =  Kx(j+2)/a2,  and  $z(t)  =  e_t  ^2 ,  the  cf  of  a 

standard  normal  random  variable.  Cearly,  if  lim^^A  (j  )  /  =0, 

then  limo^oo4>Y(t  |G)  =  <J>z(t),  or  the  limiting  distribution  of  Y 

is  standard  normal.  This  is  the  case  for  each  of  the  four 
distributions  we  have  investigated. 


r 
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Expanding  the  exponential  of  (2. A)  in  a  Taylor  series, 
we  have : 


$v(t|0)  =  <j>7(t) 


|1  +  l 

k=l 


l  A(j)(it)j+2/(oj(j+2)l) 

j  =  l 


)  k/kl|; 


(2.5) 


which  is  valid  for  jt]  <1.  The  fact  that  the  series  (2.5) 
is  not  valid  for  all  t  was  the  source  of  some  confusion  re¬ 
garding  the  convergence  properties  of  the  Edgeworth  series 
(see  Fisher  [5]).  Cramer  [2]  ended  the  claims  that  the  Edge- 
worth  series  is  convergent. 


Now,  if  we  expand  (2.5),  collecting  terms  with  the  same 

power  of  a,  we  have  a  polynomial  in  (o  ^)  whose  coefficients 
are  polynomials  in  powers  of  (it)  : 


<f>Y(t  |  0) 


1 


+  l  B.(it)a'J 

j-1  J 


where  B^it)  =  (it) 3/ 3 1  , 


(2.6) 


B2(it)  -  A2(it)4/4!  +  A2  (it)6/  (  ( 3 ! ) 2 2 ! )  ,  and 
B3(it)  -  A3(it)5/5!  +  A1A2(it)7(3!A!)  +  A2(it)9/ (3 ! )4 , 


etc . 


With  a  minor  modification  of  the  notation  used  by  Hill  and 
Davis  [7],  we  can  develop  a  general  expression  for  (it)  as 

follows . 


Denote  by  IK  ,  a  partition  of  the  po 
into  £  positive  integers: 


Oi-U-VC 


A-  V  4  A- 


l  P  •  s  .  ,  £ 

l  i 


i=l 


3  , 

(2.7) 


and  define  three  functions  of  the  partition, 

k 


m(ITj) 

a(V 


l  p.(s.+2)  =  2 £H- j  , 
i=l  1  1 

k  p . 

n  p.!((s.+2)!)  i 

i=l  1  1 


and 


(2.8) 

(2.9) 
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A(H.)  -  n  o  .  (2.10) 

J  i=l  Si 

Combining  equations  (2.8),  (2.9),  and  (2.10),  we  have  : 

B  (it)  =  l  a(n.)A(lI  )(it)2£+j;  (2.11) 

J  n.  J  J 

j 

where  the  summation  is  taken  over  all  distinct  partitions  of 
the  integer,  j . 

The  next  step  is  to  invert  the  cf  (2.6)  to  cover  the  pdf 
of  the  random  variable,  Y.  The  inversion  formula  is: 

fY(y|0)  =  (2tt)  -1/”coexp{-ity}ij>Y(t  |  0)dt .  (2>L2) 

Applying  (2.12)  to  (2.6),  which  is  valid  only  for  ]t|  <  1, 
produces  the  Edgeworth  series  for  the  pdf.  Note  that  the 
result  does  not  equal  fy(y|0)  because  the  range  of  integra¬ 
tion  exceeds  the  radius  of  convergence  of  the  expansion. 

To  apply  (2.12)  to  the  expansion  (2.6)  term  by  term,  the 
following  two  results  are  required: 

a)  From  Fourier  analysis, 

(-l)j~Vj)(y)  =  (2TT)"%/“co(it)j'1exp(-ity-t2)/2]dt;  (2.13) 

where  4>^(y)  =  the  j  th  derivative  of  a  standard  normal  cdf. 
(Note:  <f>^(y)  =  fz(y).) 

b)  From  Rodrigues'  Formula  for  Hermite  Polynomials, 

<j>  (1)He  (y  |  j  -1)  =  (-l)j~Vj)(y).  (2.14) 

Applying  (2.13)  and  (2.14)  to  (2.6),  the  powers  of  (it)  in 
(2.11)  can  be  immediately  replaced  by  the  corresponding 
Hermite  polynomial,  yielding  the  Edgeworth  asymptotic  expan¬ 
sion  for  the  pdf  of  the  random  variable,  Y: 

fY(y|0)  =  fY(y|0  n)  =■  <J>^  (y)  (1  +  l  C.  (y)o"j};  where  r':  15) 

j  =  l  J 

C,(y)  =  [  a(n.) A(n.)He(y|2£+j) .  (2.16) 

J  n_,  J  J 
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Finally,  to  get  the  cdf  ^Fy(y|G)  =  f^fy  (u  1 0)duj  ,  the  inte¬ 
gration  is  performed  term  wise  to  yield: 

FY(y|0)  *  Fy(y|o,n)  =  4>(0)(y)  -  4>(1)(y)  |  (y)o"j  (2.17) 

where  (y)  are  the  same  as  (y)  except  each  Hermite  poly¬ 

nomial  is  reduced  in  order  by  1 ,  i.e., 

D.(y)  =  l  a(n.)X(n.)He(y| 2£+j -1) .  (2.18) 

J  Hj  J  J 

3 .  Determining  a  J-fraction  Corresponding  to  the  Edgeworth 
Series 


To  generalize  the  Edgeworth  series,  let  f(y|G,o)  be  a 

function  of  y,  depending  on  a  possible  vector  parameter,  Q 
and  o.  Then  the  Edgeworth  expansion  of  the  function  is 

CO 

f (y[G,o)  =  l  P.(y|Q)o  with  P.(y|0)  polynomials.  Wall 
j=0  J  J 

(pp.  362-3)  [10],  gives  several  criteria  for  establishing  the 
existence  of  a  continued  fraction,  J-fraction  in  Wall's 
terminology,  corresponding  to  a  given  power  series.  The 
J-fraction  is  of  the  form: 


b^  +  o  -  a^ 

\)2  +  a  -  7 


If  the  Edgeworth  series  had  been  derived  directly  from  the 
Fourier-Stieltjes  transform,  without  going  through  the 
various  transformations  and  rearrangements,  the  existence  and 
convergence  of  the  J-fraction  could  be  established  in  a 
simpler  manner  (Henrici  [6]).  However,  for  the  Edgeworth 
series,  the  path  of  least  resistance  seems  to  be  to  use  a 

?uotient-difference  type  of  algorithm  (Henrici  [6]  ,  Wall 
10]),  to  determine  if  a  tractable  law  for  the  J-fraction 
coefficients  exists.  This  is  where  the  powerful  functions  of 
MACSYMA  can  be  used;  i.e.,  the  algorithms  can  be  implemented 
symbolically  and  hopefully  the  resulting  expressions  will  be 
reduceable  to  the  point  where  some  law  of  formation  can  be 
discovered.  [Note:  in  the  quotient-difference  type  algo¬ 
rithms,  the  coefficients  of  the  J-fraction  are  obtained  in 
terms  of  the  coefficients  of  the  Edgeworth  series.] 


A  starting  point  would  be  to  try  the  above  approach  on 
one  of  the  simpler  Edgeworth  expansions,  such  as  that  for  the 
incomplete  gamma  function.  Assuming  success  here  (or  fail¬ 
ure)  ,  one  could  approach  the  Edgeworth  series  for  central  and 
noncentral  beta  distributions.  If  the  law  for  the  J-fraction 
coefficients  can  be  determined  through  MACSYMA,  then  the  next 
step  is  to  try  to  establish  convergence  of  the  continued 
fraction  to  the  desired  function.  Doing  that,  the  Edgeworth 
asymptotic  expansion  can  then  be  bounded  through  the  J-frac¬ 
tion. 

4 .  Conclusions 

Perhaps  the  application  of  MACSYMA  to  the  problem  of 
bounding  Edgeworth  asymptotic  expansions  will  not  produce  any 
good  results.  In  any  case,  it  can  be  no  less  successful  than 
the  many  previous  attempts  by  other  methods.  It  is  hoped 
that  the  occasional  orderliness  of  statistical  distribution 
functions  will  surface  in  this  attempt. 
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ABSTRACT 

Statisticians  have  recently  been  concerned  with  robust  regression  concepts. 
Introduction  of  M-estimation  procedures  somewhat  extends  the  computational 
burden  over  ordinary  multiple  linear  models  with  normal  errors,  but  does  not  yet 
begin  to  test  the  computational  capabilities  of  modern  computer  architectures. 
Consider,  for  contrast,  a  reasonably  straightforward  application  in  image 
processing.  An  image  may  be  represented  as  a  function,  f:  R^->-R,  taking  the 
real  plane  into  the  real  line,  z  «*  f(x,y)  being  the  intensity  at  point  (x,y) 
and  x  and  y  be  respectively  the  horizontal  and  vertical  locations  of  a  point 
in  the  image.  A  very  simple  noisy  picture  may  be  represented  as  a  regression 
problem  by 

Zij  "  f(xi*  V  +  Eij  (1) 

where  the  are  the  usual  white  Gaussian  noise.  This  problem  corresponds 

to  a  snowy  TV  picture  in  weak  reception  area.  Since  f  in  general  is  a  strongly 
nonlinear  function  of  x^  and  y^  it  is  clear  that  some  nonlinear  nonparar.x-t: ic 

regression  methodology  is  necessary  even  in  this  simple  case.  Perhaps  the 
do ■  st  work  to  solving  this  problem  is  that  of  Wahba  (1979) .  See  also 
Wegman  and  Wright  (1980) .  Even  this  spline  approach  is  not  entirely  adequate 
since  splines  are  required  to  be  smooth  and  an  image  may  have  very  sharp 
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discontinuities  (high  contrast).  In  certain  satellite  remote  sensing  appli¬ 
cations,  imaging  sensors  operate  at  very  low  light  levels  where  individual 
photons  are  counted.  The  white  Gaussian  noise  assumption  is  then  replaced 
by  a  Poisson  noise  distribution.  Moreover  since  these  light  levels  are  in 
a  nonlinear  response  region  for  film  or  other  imagers  (due  to  so  called 
reciprocity  failure),  the  additive  noise  assumption  is  no  longer  appropriate. 
Thus  a  regression  model  might  take  the  more  complex  form 


zij  “  f(xi*  yi>  0  eij  (2) 

where  o  is  a  nonlinear  binary  operator  and  (to  coin  a  phrase)  is  dark 

non-Gaussian  noise.  To  complicate  the  picture  even  more,  we  could  ask  for 
color  images  which  means  that  f  must  be  a  vector  valued  map,  e.g., 
f:  r2-*R-*.  We  might  model  this  as 


where 

image. 


(rij>  g^)  *=  f(x±>  Yj)  o  e±i 


(3) 


rij ’  bij  an<*  rePresent  the  red,  blue  and  green  components  of  the 
Notice  that  r ,  b^  and  g^  will  in  general  be  correlated  random 


variables.  One  further  complication  is  to  suggest  we  might  be  interested  Jn 
a  motion  picture.  We  thus  introduce  a  time  series  aspect  to  our  evolving 
model  which  now  may  look  like 


(rijk’  bijk*  8ijk)  “  f(xi’  yi*  tk)  0  Eijk  * 
and  clearly  f:  R^+R^.  To  put  some  simple  numbers  to  this  example  will  help 
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clarify  our  point.  Even  discounting  the  computational  burden  required  to 
estimate  f  in  the  first  place,  a  high  resolution  image  normally  is 
represented  digitally  on  1024  x  1024  grid  (of  pixels) .  Thus  each  image 
requires  2^®  evaluations  of  f  (itself  a  vector-valued  function) .  For  a 
simple  20  minute  film  there  are  28,800  images.  If  we  were  doing  the  digital 
processing  image  by  image,  this  would  mean  we  had  28,800  nonlinear,  non- 
parametric  non-Gaussian  regressions  to  do  and  to  evaluate  each  of  the  28,800 
vector-valued  function  estimators  at  2^0  points.  All  this  is  required  to 
process  digitally  20  minutes  of  color  film  taken  in  poor  light,  a  relatively 
realistic  assignment.  The  dimensionality  and  sample  size  requirements 
clearly  demonstrate  the  need  for  innovative  development  of  statistical 
algorithms  based  on  a  sound  knowledge  of  modern  computer  architectures. 

Some  recent  developments  in  microelectronic  technology  have  revolution¬ 
ized  computer  design.  Very  large  scale  integrated  circuit  technology  (VLSI)  has 
revolutionized  the  concept  of  central  processing  units.  VLSI  circuit  chips 
now  can  contain  a  multiplier  which  makes  parallel  and  network  arrangements  of 
processors  possible  in  a  relatively  inexpensive  fashion.  Processors  may  be 
connected,  for  example,  not  only  in  parallel  arrays  but  in  orthogonal  or 
hexagonally  connected  arrays.  Such  innovative  computer  architectures  allow 
for  a  totally  different  approach  to  algorithm  development.  Some  examples 
using  matrix  manipulations  suggest  that  statistical  algorithms  for  multi¬ 
dimensional  data  can  be  formulated  in  a  fundamentally  different  way.  See, for 
example.  Rung  and  Leiserson  (1978)  or  Mead  and  Conway  (1980).  Ve  advocate 


an  organic  approach  to  algorithm  development.  That  is,  rather  than  having 
a  theoretician  develop  a  formula  which  is  then  translated  by  a  programmer 
to  a  computer  algorithm,  we  think  there  is  much  to  be  derived  from  an 
Integrated  approach. 
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ABSTRACT 

In  the  study  of  probabilities  of  large  samples  one 
often  encounters  integrals  of  the  form 

b 

(1)  I  =  /  <(>  (x)  [  f  (x)  ]ndx, 

n 

a 

where  <j>(x)  and  f(x)  are  continuous  functions  defined  on  the 
finite  or  infinite  interval  [a,b]  and  f(x)  is  positive  there. 
Long  ago,  Laplace  made  the  observation  that  the  major  con¬ 
tribution  to  the  integral  should  come  from  the  neighbourhoods 
of  the  points  where  f(x)  attains  its  greatest  value. 
Furthermore,  he  showed  that  if  f(x)  attains  its  maximum 
value  only  at  the  point  £  in  (a,b)  where  f'  (£)  =  0  and 

f '  '  (£  )  <  0,  then  as  n  -*•  “ 

(2)  I„  t^Ter)''- 

This  formula  is  now  known  as  the  Laplace  approximation. 


However,  it  is  not  infrequent  to  come  across  integrals 


of  the  form  I  that  do  not  satisfy  the  conditions  necessary 
for  the  validity  of  the  approximation  in  (2).  Hence 
extensions  of  and  modifications  to  the  method  of  Laplace 
must  be  made  in  order  to  obtain  the  behaviour  of  these 
integrals  for  large  values  of  n.  Here  we  shall  consider 
two  such  cases.  Each  of  these  is  illustrated  by  a  specific 
example . 

The  first  modification  concerns  the  integral 

2 


(3) 


I  (n) 


I 


-  X 

xe 


r  1  +  6  (x)  -,n 
2 


dx , 


--  CO 


where  fl(x)  is  given  by 

it  __  2 

( 4  )  0  ( x )  =  — £  /  «.  U  d  u . 

/it  0 

This  integral  occurred  in  a  problem  in  probability  theory. 
Note  that  the  function  *5fl  +  0(x).]  mono  t  or.  i  ca  1  1  y  increases 
from  0  to  1  as  x  varies  from  -<*>  to  +°°.  Hence  the  greatest 
value  of  this  function  is  not  attained  at  a  fir i te  point 
but  at  infinity,  and  the  conditions  for  the  La  place  approxi¬ 
mation  are  violated.  Nevertheless,  we  shall  snow  that  J(n) 
has  the  behaviour 


(5)  1 ( n ) 


/tt  log  (n-t-i  )_ 
(nM)  '  ~ 


/ 1  log  log  ( n  +  1 )  _ 1 _ . 

4  .  - : — TTT  (  n  1 )  log  (n+1  ) 

(n+l)/logvn+l) 


as  n 


The  second  modification  deals  with  the  incomplete 
Beta-type  integral 


i 


the  integral  (6)  becomes 

1 

(11)  J(n)=  j  <f>(t;n)if(t)]n  dt, 

a 

n 

which  is  indeed  of  the  form  in  (1) .  Note  that  the  function 
<(>  and  the  lower  limit  a  now  depend  on  n  and  hence  the  Laplace 
approximation  (2)  does  not  apply  directly.  However,  we  shall 
use  a  modification  of  Laplace's  method  to  show  that  J(n) 
has  the  asymptotic  expansion 

(12)  J(r>  [  f  ( a )  ] n  {  -1  -'-n  +  -2-~n  +  as  n  -*•  <*> , 

n  2 

n 

where  the  coefficients  c,  are  bounded  functions  of  n.  As 

i ,  n 

an  application  of  this  result,  we  show  that  the  tail  prohr - 


bility  of  the  sample  p-quantile  decays  exponentially. 
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List  of  Lectures 


DATE 

SPEAKER 

TOPIC 

Oct . 

7,  1980 

J.  Harrison 

Flows  without  periodic  orbits 

Oct . 

24,  1980 

H.-O.  Peitgen 

Continuation  and  bifurcation 

Feb . 

26,  1981 

B.  C.  Eaves 

An  introduction  to  solving  equations 
with  P.L.  homotopies 

Mar . 

24,  1981 

H.  Keller 

Continuation  methods  in  seismic 
ray  tracing 

Mar. 

25,  1981 

S.  Smale 

Fundamental  theory  of  algebra: 
historical  remarks  and  new 
perspectives  from  computer  science 

Mar. 

26,  1981 

S.  Smale 

Cost  of  finding  zeros  of  maps  by  a 
variation  of  Newton's  method 

Mar. 

27,  1981 

H.  Keller 

The  optimization  problem  in 
continuation  and  homotopy 

Mar. 

31,  1981 

L.  Watson 

Homotopy  methods — the  Colt  45  of 

nonlinear  equation  solvers 


